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Finite  Difference  Computation  of  Wave  Drag  and  Pressure  on 
Slender  Bodies  of  Revolution  at  Transonic  Speeds 

with  Zero-Lift* 


Li  Hsiu-ying  and  Lo  Shih-chiin 
(Northwestern  Polytechnical  University) 


ABSTRACT 

A  transonic  axisymmetrical  potential  equation  with  large  dis¬ 
turbance  in  the  free  stream  direction  and  small  disturbance  in  the 
transverse  direction  is  solved  by  using  the  Murman-Cole  schemes  of 
finite  differences.  The  boundary  condition  on  the  body  is  trans¬ 
ferred  to  body  axis.  The  boundary  condition  at  far  field  is  approxi¬ 
mated  by  that  at  infinity. 

Finite  difference  equations  for.  the  potential  are  solved  by 
line-overrelaxation  along  the  radius  with  Seidel  iteration. 

In  order  to  calculate  the  pressure  on  the  body  surface,  the  potential 
is  interpolated  by  the  slender-body  theory.  The  pressure  coefficient 
is  calculated  by  the  exact  Bernoulli's  equation.  The  zero-lift  wave- 
drag  coefficients  are  obtained  by  integrating  the  pressure  coefficients 
on  the  body  surface. 

The  computational  results  for  seven  different  configurations  agree 
well  with  the  known  wind  tunnel  test  data  as  shown  in  Figure  2  to 
Figure  6.  The  experiences  obtained  from  investigation  of  mesh  spacings, 
initial-fields,  iterative  methods,  relaxation  factors,  etc.,  in 
relation  to  the  convergence  and  rate  of  convergence  may  be  interesting 
to  engineers. 

A  linearized  analysis  of  the  stability  and  the  convergence  in  line 
overrelaxation  of  difference  equations  for  steady  axisymmetric  small 
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perturbation  potential  flow  is  carried  out  and  the  conclusions  are 
shown  in  Table  1.  The  numerical  computations  do  agree  with  the 
theoretical  conclusions. 

It  must  be  pointed  out  that  as  M»  is  very  close  to  unity  or  the 
mesh  spacing  is  shortened,  the  iterative  computation  does  not  con¬ 
verge  to  the  usual  degree  of  accuracy.  This  fact  might  be  explained 
by  a  decrease  in  the  artificial  viscosity  in  the  potential  difference 
equation  at  locally  supersonic  points. 

The  mesh  is  taken  to  be  62  x  16,  and  the  initial  field  for  the 
iteration,  0.  The  numbers  of  iterations  necessary  before  conver¬ 
gence  is  achieved  are  approximately  150,  40  and  300,  corresponding  to 
M-<  1 .  M->  1  and  1  ,  respectively.  The  relaxation  factors 

are  taken  to  be  0.9<a>»<1.7,  O.S^cb^i.Oj  for  M-<C  1 

and  0.8<£a>»<0.9,  0.8<«»<0.9,  for  M->  1  .  Here 

are  the  relaxation  factors  of  locally  subsonic  and  locally  supersonic 
points,  respectively. 

I .  INTRODUCTION 

According  to  the  area  law  for  transonic  speeds,  the  zero-lift 
wave  drag  of  an  aircraft  flying  at  transonic  speeds  is  equal  to  that 
of  its  equivalent  body  of  revolution  under  given  conditions.  Therefore, 
it  is  of  practical  value  to  aircraft  design  and  analysis  to  present 
a  convenient  method  and  procedure  for  calculating  the  wave  drag  of  bodl 
of  revolution  at  transonic  speeds  with  zero-lift. 

The  mixed  finite  difference  scheme  ^  is  an  effective  method 
for  solving  transonic  constant  steady  potential  flow  problems  with 
computers.  In  References  [2],  [3]  and  [4],  the  canonical  small  per¬ 
turbation  velocity  potential  equation  and  boundary  conditions  have 
been  used  to  satisfactorily  compute  the  pressure  distribution  about 


a  moving  parabolic  body  of  revolution.  However,  there  is  as  yet  no 
publicated  material  on  the  computation  of  wave  drag  and  pressure 
distribution  about  a  body  of  revolution  of  any  given  shape.  Keep¬ 
ing  in  mind  the  requirement  of  practical  applicability,  we  have  used  the 
mixed  finite  difference  scheme  to  compute  and  study  the  transonic 
axisymmetric  potential  flow  about  a  body  of  revolution  of  any  shape. 


The  canonical  small  perturbation  velocity  potential  equation  has 
a  singularity  at  the  nodal  point  on  the  bow  of  an  object  with  a  blunt 
bow,  and  is  not  applicable  at  other  nodal  points  either.  Hence,  it 
cannot  be  applied  to  an  object  that  may  take  on  any  shape.  In  Reference 
[5] ,  a  potential  equation  with  large  perturbation  in  the  free 
stream  direction  and  small  perturbation  in  the  transverse  direction 
and  the  corresponding  boundary  conditions  have  been  adopted.  This 
avoids  the  singularity  of  the  canonical  small  perturbation  equation  at 
the  nodal  point  on  a  blunt  bow.  In  theory,  this  can  increase  the 
accuracy  of  the  computation  of  the  flow  in  the  vicinity  of  the  nodal 
point.  Furthermore,  it  does  not  require  any  additional  computing  time. 
Based  on  the  above  considerations,  we  also  adopt  this  equation  and  its 
boundary  conditions  in  this  work. 

The  finite  difference  equations  are  solved  by  line-overrelaxation 
along  the  radius  with  Seidel  iteration.  Stability  and  convergence  are 
important  criteria  of  a  successful  iterative  computation.  To  provide 
a  guideline  for  the  iterative  computation,  we  have  performed  an  analysis 
of  stability  and  convergence.  Owing  to  the  complexity  of  the  nonlinear 
problem,  this  analysis  has  been  carried  out  under  the  condition  of 
local  linearization. 

II.  FUNDAMENTAL  EQUATIONS  AND  FORMULAS 

Take  the  Ox-axis  to  be  along  the  axial  direction  of  the  object, 
and  the  Or-axis  to  be  along  the  radial  direction  of  the  object.  The 
origin  of  the  coordinate  system  is  taken  to  be  at  the  apex  of  the  bow 
(Figure  1) .  Assume  that  the  large  perturbation  is  along  the  x-direction 
while  the  small  perturbation  is  along  the  r-direction.  The  perturbing 
velocity  potential  equation  is  given  by 


The  boundary  condition  at  far  field  is  approximated  by  that  at  infinity 


«P  —  0 


(5) 


Along  the  x-axis  in  front  of  the  object,  the  flow  is  symmetrical  with 
respect  to  the  x-axis. 


cp, -<?>.-  0 


(6) 


Let  Ar ^  -  2Ar^.  Inserting  the  boundary  condition  on  the  surface 
of  the  object  into  the  velocity  potential  equation,  one  obtains  the 
expression  (Refer  to  Figure  1) : 

~  Wi.x-v-.x-RR'  Cg.+«JP,(  x ,  rjj>  +  o(aH) 


Here,  j  and  k  are  the  index  numbers  of  the  mesh  along  the  x-direction 
and  r -direction,  respectively. 

Inserting  the  boundary  condition  given  by  Equation  (4)  for  the 
apex  of  the  bow  into  the  velocity  potential  equation,  one  obtains  the 
expression 


Here,  jl  is  the  j-index  of  the  apex  of  the  bow. 

Inserting  the  boundary  condition  on  the  x-axis  in  front  of  the 
object  into  the  velocity  potential  equation,  one  obtains  the  expression 


At(*m-*w>+  0  (v) 


Use  the  mixed  finite  difference  scheme,  and  substitute  all  the 
derivatives  in  Equation  (1)  and  the  equations  obtained  by  inserting  the 
boundary  conditions  in  Equation  (1)  with  the  corresponding  finite 
difference  expressions.  One  then  obtains  a  set  of  finite  difference 
equations.  To  linearize  this  set  of  equations,  regard  the  coefficients 
of  all  the  terms  as  known  quantities,  and  use  an  iterative  method 
of  computation.  The  finite  difference  equations  are  solved  by  line- 
overrelaxation  along  the  radius  with  Seidel  iteration.  The  overrelaxati 
operation  is  defined  by 

(7) 


in  which  is  the  result  of  the  nth  iteration  without  over- 

3  i K 

relaxation  treatment,  and  ta  is  the  relaxation  factor. 

The  pressure  coefficients  on  the  surface  of  the  object  are  given 
by  the  exact  Bernoulli's  equation: 


where  <|>x  is  derived  from  the  slender  body  theory  ; 


(H*  *+ RRr)gmln  (R/rt) 


The  pressure  coefficients  at  r=constant  in  the  flow  field  are 

r  2  ff .  /  Y  i  •)  w»  (  Op,)?*  .  (g*.)*,*,  ■  V  T  -  i  -  I  ) 

_(YUM-V  9-  +  2<7i  2  ql  f  (9) 

The  pressure  drag  coefficient  (not  including  bottom  drag) 
is  given  by  the  following  expression. 

^  _  2*  f2*  n  uof  /ia\ 


'-“-T-J.  C’RR'‘ 


In  the  above  equation,  s  is  the  reference  area,  and  is  the  object 
length  cut  off  by  the  tail  branch  rod. 


The  zero-lift  wave  drag  is  defined  as 


“*C«» 


where  M*  is  the  critical  Mach  number. 


III.  COMPUTED  EXAMPLE 


The  examples  computed  consist  of  seven  slender  bodies  of 
revolution  with  different  configurations.  The  surface  pressure  distribu 
tion  and  zero-lift  pressure  drag  coefficients  of  a  parabolic-arc  body  of 
revolution  having  fineness  ratio  of  12  are  given  in  Figures  2(a) 
and  2(b) ,  respectively.  Figure  3  gives  the  zero-lift  wave  drag 
coefficients  for  a  sing-body  combination  and  those  for  the  body  alone. 
The  body  consists  of  a  cylinder  and  a  missle  ty£~  forebody,  and  the  wing 
is  a  delta  wing.  Figure  4(a)  gives  the  zero-lift  wave  drag  coeffi¬ 
cients  for  a  swept-wing  and  curved-body  combination^  and  those  for 
the  body  alone.  Figure  4(b)  gives  the  computed  pressure  distribution  of 
a  body  of  revolution  equivalent  to  the  swept-wing  and  curved-body  com¬ 
bination.  The  shock  wave  positions  obtained  from  Figure  4(b)  are 
given  in  Figure  4(c).  Figures  5(a)  and  6(a)  give  the  cross-sectional 
area  distribution  in  the  axial  direction  for  a  certain  aircraft  I  and 
aircraft  II,  respectively.  The  zero-lift  wave  drag  coefficients 
corresponding  to  these  area  distributions  are  shown  in  Figures  5(b) 
and  6(b),  respectively.  Except  for  the  results  shown  in  Figure  4(b), 
all  the  computed  results  have  been  compared  with  the  wind  tunnel  test 
data.  There  is  a  fairly  good  agreement  between  the  computed  and 
experimental  results.  In  Figure  5(b),  the  results  obtained  from 
finite  difference  computations  are  also  compared  with  those  obtained 
from  the  linearized  theory.  It  is  obvious  that  the  former  is  superior 
to  the  latter. 
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Figure  2(a)  Pressure  distribution  on  surface  of  a  parabolic-  ,g. 

arc  body  of  revolution  having  fineness  ratio  of  12 ' 
Keys  1-  Finite  difference  computation;  2-  Experimental  results 
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Figure  2(b)  Zero-lift  pressure  drag  coefficients  of  a  parabolic-arc 
body  of  revolution  having  fineness  ratio  of  12(6) 

Key:  1-  Finite  difference  computation;  2-  Experimental  results 


The  discrepancy  between  the  results  computed  by  finite  difference 

method  and  those  obtained  experimentally  is  chiefly  due  to  two  factors. 

First,  in  the  finite  difference  computation,  the  effect  of  the  boundary 

layer  has  not  been  taken  into  consideration.  Secondly,  the  experimental 

data  are  subject  to  errors  due  to  their  dependence  on  the  extent  the 

r  o  1 

wind  tunnel  is  choked1  1  . 

The  values  of  M°°  in  the  examples  range  from  0.8  to  1.4.  There 
are  from  11  to  13  Moo  for  each  example.  Computer  time  required  on 
the  Felix  C-256  machine  is  around  two  hours  while  it  is  about  a  half 
hour  on  the  Siemens  7760  machine. 
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Figure  3.  Zero-lift  wave  drag  coefficients  for  a  body  of  revolution 
equivalent  to  a  delta-wing-body  combination  and 

for  the  isolated  body  consisting  of  a  cylinder  and  a 
missile  type  forebody. 

Key:  1-  Finite  difference  computation;  2-  Experimental  result; 

3-  Wing-body  combination;  4-  Isolated  body. 
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Figure  4(a).  Zero- lift  wave  drag  coefficients  for  a  body  of  revolution 
equivalent  to  a  swept-wing  and  curved-body  combination  l'] 
and  for  the  isolated  body. 

Key:  1-  Finite  difference  computation;  2-  Experimental  result; 

3-  Wing-body  combination;  4-  Isolated  body. 
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Figure  4(b). 


Computed  pressure  distribution  of  a  body  of  revolution 
equivalent  to  a  swept-wing  and  curved-body  combination 
17]  at  M_  =  0.98. 
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Figure  4(c).  Shock-wave  positions  obtained  from  Figure  4(b). 
Key:  1-  Finite  difference  computation;  2-  Experimental  result 


Figure  5(a).  Cross-sectional  area  distribution  of  the  aircraft 


Figure  5(b).  Zero-lift  wave  drag  coefficients  of  the  aircraft. 

Key:  1-  Finite  difference  computation;  2-  Experimental  result; 
5-  Linearized  theory. 
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Figure  6(a).  Cross-sectional  area  distribution  of  the  aircraft  I 
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Figure  6(b)  Zero-lift  wave  drag  coefficients  of  the  aircraft  II. 
Key:  1-  Finite  difference  computation;  2-  Experimental  result 


IV.  DISCUSSION  OF  METHODS 

Comparison  of  computed  results  with  experimental  data  shows  that 
it  is  better  to  treat  the  apex  of  the  bow  and  the  end  point  of  the 
tail  as  points  on  the  inner  side  of  the  surface. 

To  study  the  effect  of  mesh  spacing  on  the  accuracy  and  conver¬ 
gence  of  the  results  of  computation,  we  compared  the  results  obtained  by 
distributing  40,  80  and  160  mesh  points,  respectively,  on  the  object 
along  the  axial  line.  We  found  that  the  choice  of  40  mesh  points 
along  the  axial  line  allows  faster  convergence,  saves  computer  time 
and  meets  the  requirement  on  accuracy  of  the  job.  Therefore,  our 
choice  of  the  mesh  is  as  follows:  The  distance  on  the  x-axis  between 
x»0,  and  x»Ld  is  divided  into  40  equal  spacings.  To  the  front  of 
x=0,  10  spacings  are  “-^ken,  with  the  first  two  each  equal  to  Lp/40,  and 
the  remaining  step  .  vs  increasing  in  such  a  manner  as  to  form  a 
geometric  progress i  a  common  ratio  of  2.  To  the  back  of 

x»Ld,  11  spacings  a  t  ,  with  the  same  variation  in  step  length 

as  that  for  the  points  ta^en  in  front  of  x»0.  Fifteen  spacings  are 
taken  along  the  r-axis,  with  Ar2*  Ar3“  2Ar^  and  the  remaining  step 
lengths  forming  a  geometric  progression  with  a  common  ratio  of  2. 
Experience  shows  that  Ar^  can  have  any  value  lying  within  the  range 
Ar^  -*-0.  He  take  it  to  be  LD/80.  Hence,  the  total  number  of  mesh 
points  is  62  x  16,  and  the  farfield  boundaries  are  x«-12.8LD,  26.6LQ 
and  r*204.8Ln. 


We  have  also  carried  out  trial  computations  for  different  far- 
field  boundaries.  The  results  show  that  the  requirement  on  accuracy  of 
the  job  is  satisfactorily  met  by  taking  them  at  10-20  in  front  of  and 
to  the  back  of  the  object  along  the  x-axis,  and  at  100  LQ  along  the 
r-direction.  However,  it  may  be  well  in  the  case  of  M,»  =1  to  take 
the  boundary  points  farther  out  along  the  r-direction.  With  this 
choice  of  the  boundary  points,  the  farfield  boundary  condition  may  be 
approximated  by  Equation  (5) .  If  the  farfield  boundaries  are  taken  to 
be  too  close  to  the  object,  then  there  will  be  larger  errors  produced  by 
approximating  the  farfield  boundary  condition  with  Equation  (5) . 

We  have  compared  the  results  obtained  by  using  zero  initial-field 
and  those  obtained  by  taking  the  initial  field  to  be  the  final  field 
computed  from  a  neighboring  Moo  value  that  has  been  previously 
obtained  from  computation.  The  results  show  that  the  former  con¬ 
verges  faster.  While  the  latter  also  converges,  it  requires  many 
more  iterations  to  attain  the  same  degree  of  accuracy.  Therefore,  the 
initial  field  is  taken  to  be  zero  in  this  work. 

Both  computational  work  and  analysis  of  the  local  linearization 
theory  (Section  V)  show  that,  under  supercritical  conditions,  the  use  of 

the  simple  iterative  method  leads  to  instability  in  the  computation. 

.  2 

Except  for  the  velocity  determinant  1-M  for  which  the  simple 
iterative  method  has  been  used,  we  use  the  improved  Seidel  iteration 
throughout  this  work. 

The  results  of  the  computations  show  that  the  optimal  relaxation 
coefficient  depends  on  M,,,  and  the  step  length.  For  the  same  mesh 
spacing,  the  closer  M«  is  to  1,  the  smaller  the  value  of  the  relaxation 
factor  should  be  taken  within  the  range  given  below.  Otherwise,  the 
computation  may  not  converge.  For  the  same  M«  ,  the  smaller  the  step 
length,  the  smaller  the  optimal  relaxation  coefficient.  This  is  in 
agreement  with  the  conclusion  obtained  from  the  stability  analysis 
(Section  V) .  In  this  paper,  the  relaxation  coefficients  for  the 
62  x  16  mesh  are  taken  from  the  following  ranges: 

For  M  <  1,  0  >7|  o.9<o,<1.0» 

for  M  >1,  0.8<ffli<0.9,  0.8<o,<0.9 


Here,  oj^  and  co^  denote,  respectively,  the  relaxation  factors  for  the 
locally  subsonic  points  and  locally  supersonic  points. 


For  M»<1,  the  convergence  is  usually  monotonous.  In  general, 
after  about  150  iterations,  AV<10*4  (A«P*  1a<p1?./a<p.11.>.I  denotes 

the  maximal  relative  iterative  error  of  the  velocity  potential  field  . 
For  Moo  >1,  the  convergence  goes  through  an  oscillatory  phase  before 
it  becomes  monotonous.  In  general,  A<1P<10'4  can  be  attained  after  40 

iterations.  For  Mco*\*l,  individual  points  may  cause  the  velocity 
potential  not  to  converge  to  the  usual  degree  of  accuracy.  However, 
in  general,  after  200  to  300  iterations,  the  aerodynamic  coefficients 
can  be  made  to  converge.  Similar  results  are  obtained  when  the  mesh 
is  shortened  along  the  x-direction.  This  may  be  due  to  a  decrease 
in  the  artificial  viscosity  in  the  finite  difference  equation  at 
locally  supersonic  points.  Said  artificial  viscosity  is  given  by 


( 1  -M*)<*>«A* 

Both  M  -*•  1  and  Ax-*-o  cause  the  artificial  viscosity  to  become  smaller, 
making  it  difficult  for  the  computation  to  converge.  This  problem 
requires  further  study. 


V. 


STABILITY  AND  CONVERGENCE  UNDER  THE  CONDITION  OF  LOCAL  LINEARI ZATIQl 


When  Equation  (1)  is  solved  with  the  iterative  method,  the  non¬ 
linear  terms  in  the  equation  have  actually  been  locally  linearized  at 
every  point. 

Assume  that  the  coefficient  1-M  of  the  <l>xxterm  in  Equation  (1) 
is  a  constant.  Then,  at  the  locally  subsonic  and  locally  supersonic 
points.  Equation  (1)  can  be  transformed  respectively  into 


±<P„+(r<l>,),/r-  0 


(ID 


Let  the  step  lengths  Ax  and  Ar  be  constants.  Then  the  finite 
difference  equations  for  Equation  (11)  to  be  solved  by  line-over¬ 
relaxation  along  the  radial  direction  becotne 


A.*1 

-Stt- 29, 
Ax* 


-  <•)  _ a*  do  i 

.  «  * 


Ar* 


fi.H 


r*2Ar 


-  0 


(12a) 

(12b) 


Equations  (12a)  and  (12b)  are  equivalent  to  the  following 
differential  equations,  respectively: 


<0 


0 


(12a' ) 
(12b') 


Equations  (12a')  and  (12b')  are,  respectively,  hyperbolic  and 
parabolic,  and  have  inherent  stability  problems.  The  solutions  of 
Equation  (12a)  and  (12b)  can  be  expressed  in  the  form 
9(x,  r,  i)+5(x,  r,  O  ,  where  <t>(x,r,t)  is  the  exact  solution  and 
6(x,r,t)  is  the  error.  5(x,r,t)  also  satisfy  Equations  (12a)  and 
(12b) ,  and  can  be  expressed  as 

b(x,  r,  t  (13) 

In  the  above  equation,  i  =  /^T,  ^and  B2  are  any  real  numbers  and 
a  is  a  function  of  and  B2*  The  condition  of  stability  requires 
that  the  error  not  increase  with  time,  i.e., 

n<i  (i4) 


Making  use  of  Equations  (7),  (12a)  and  (13),  we  obtain 
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(15) 


The  condition  under  which  the  above  equation  satisfies  the  inequality 
(14)  is  0<®<a>,<2  .  Here,  decreases  with  decreasing  step 
length . 

Convergence  refers  to  the  convergence  of  the  solutions  for  the 
artificial- time-dependent  differential  equations  (12a')  and  (12b') 
to  the  solutions  for  Equation  (11) . 

Let  , +-i - .  Equation  (12a')  can  be  written  as 

2  A*  * 


The  solution  for  Equation  (16)  can  be  found  by  the  method  of 
separation  of  variables 

«P(*,  r,  J  (x,  r)  (17) 

m  -  l 

Here  Gm(x,r)  is  the  eigenfunction  of  the  boundary  value  problem 
equation  below. 

<r.+C„+ (?,/r+xl<7-  0 


is  the  corresponding  eigenvalue  with  A*<Aj<—<Ai.,<Ai<— t  <?#(*, 

is  the  solution  for  Equation  (11).  A  and  B  are  constants. 

m  m 
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Obviously,  when  0  <  (u  <2, 


limV(x,  r,  t  r) 

All  the  results  similarly  derived  are  listed  in  Table  1.  Compu¬ 
tational  study  shows  that  these  are  the  necessary  conditions  for 
stability  and  convergence. 


TABLE  1.  Stability  and  Convergence  Conditions  for  r-line  Over¬ 
relaxation  Iterations 


m  *  *  « 


4- a  a  <* 
S  *  a  « 


0  <!»<<»,•<  2  •  <»!<»<  I  0  <®<o>,#<  2  /.  ft7B2 


•  <®<J 


•  <®<1 


7.  • 


Keys  1-  Improved  iteration;  2-  Simple  iteration;  3-  Condition; 
4-  Stability;  5-  Convergence;  6-  Always  unstable; 

7-  oj^  decreases  with  decreasing  mesh  step  length 


VI.  CONCLUSION 

We  have  given  in  this  paper  a  method  for  computing  the  zero-lift 
wave  drag  and  pressure  distribution  on  slender  bodies  of  revolution  at 
transonic  speeds.  The  results  of  the  computations  agree  fairly  well 
with  experimental  results.  In  practice,  it  is  only  necessary  to 
take  40  mesh  points  on  the  object  along  its  body  axis.  When  the  bounda^ 


v'-JO.-V 


condition  for  the  farfield  is  taken  to  be  zero,  the  boundaries  of  the 
farfield  may  be  taken  to  be  at  10-20  times  the  length  of  the  object 
along  the  x-direction  on  either  end  of  the  object,  and  100  times 
the  object  length  along  the  r-direction.  The  finite  difference 
equations  are  solved  by  line  over-relaxation  iterations  along  the 
r-direction.  To  ensure  stability  and  convergence  of  the  iterative 
computation,  the  relaxation  coefficients  for  locally  subsonic  points 
are  taken  to  be  0<a»<2  .  For  small  step  lengths  or  when  M<»  is  close 

to  1,  the  optimal  value  of  oi  is  small.  The  relaxation  coefficients  for 
locally  supersonic  points  are  taken  to  be  a><i  and  close  to  1.  When  M» 
approaches  1,  the  iterative  computation  does  not  -on verge  to  the  usual 
degree  of  accuracy  at  individual  points.  In  this  case,  the 
standards  for  convergence  of  the  aerodynamic  coefficients  may  be 
used. 

The  problem  associated  with  difficulty  in  convergence  at 
individual  points  of  the  velocity  potential  when  Moo-*  1  requires 
further  study. 
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Characteristics  of  Metal  Fatigue  Propagation  in  the 
Initial  Stage  and  Its  Influencing  Factors 

Beijing  Aeronautical  Material  Institute  Yan  Minggao 

ABSTRACT 

In  this  paper,  the  mechanism  and  characteristics  of  fatigue 
crack  propagation  of  metals  and  alloys  near  the  threshold  region 
were  reviewed. 

According  to  the  results  of  microfractographic  analysis  of 
»  a  series  of  specimens  and  structures ,  the  predominant  fracture 

morphology  is  a  type  of  facet  or  river  pattern  in  the  near¬ 
threshold  region.  The  correlation  between  the  direction  of  the 
facet  of  different  metals  and  alloys,  and  the  crystal  lattice 
structure  could  be  described  by  the  slip  mode  and  dislocation. 

In  this  paper,  the  relation  between  fatigue  limits  and 
threshold  values  of  plate  and  notched  specimens,  as  well  as  the 
fatigue  crack  propagation  mechanism  of  a  short  crack,  were 
described  and  commented.  The  crack  propagation  behavior  of  a 
short  crack  can  in  general  be  expressed  as: 

A/C,*=*  /  (  a  )  (A ) 

where  f (a)  is  a  function  of  the  crack  length  and  specimen 
geometry . 

In  this  paper,  the  effect  of  stress  ratio,  microstructure, 
and  environment  on  the  fatigue  crack  propagation  behavior  in  the 
near-threshold  range  was  reviewed.  The  concept  of  "oxide  in¬ 
duced"  and  "roughness  induced"  crack  closure  effect  was  presented 
It  may  be  used  to  explain  the  effect  of  the  aforementioned  fac¬ 
tors  on  the  fatigue  crack  propagation  behavior  at  the  near¬ 
threshold  range  in  structural  materials. 

I .  INTRODUCTION 

Since  the  incidents  of  fractured  truck  axles  in  England 
early  in  this  century,  the  study  of  fatigue  fracture  has  at¬ 
tracted  a  wide  range  of  attention.  Furthermore,  a  great  deal  of 
experimental  research  and  theoretical  analysis  had  been  conducted 
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Nevertheless,  serious  incidents  still  occurred  from  time  to  time; 
from  the  rupture  of  the  ship  "Freedom"  in  the  forties  to  the  dis¬ 
integration  of  the  airliner  "Comet";  from  the  breaking  up  of  a 
rocket  shell  to  the  collapse  of  the  oil  drilling  rig  on  the  North 
Sea  in  the  eighties.  This  shows  that  the  inception  and  propaga¬ 
tion  of  fatigue  cracks  at  places  where  stress  is  concentrated  on 
a  structural  piece  will  be  unavoidable.  Therefore,  the  study  of 
the  initial  propagation  of  fatigue  cracks  using  notched  and 
microcracked  structural  pieces  has  important  theoretical  and 
practical  meaning  in  the  prevention  of  future  disasters. 

It  is  generally  believed  that  the  propagation  of  fatigue 
crack  can  be  divided  into  three  stages.  During  the  initial 
propagation  period,  it  is  propagating  at  an  extremely  low  rate 
(da/dN< 10_5mm/week) .  The  region  is  the  discontinuous  propaga¬ 
tion  zone.  A  type  of  crystallographic  facet  or  river  pattern 
fracture  begins  to  emerge.  When  da/dN>10-5mm/week,  then  it 
reaches  a  continuous  propagation  region.  The  fracture  morphology 
is  predominantly  fatigue  patterns.  When  da/dN> 10~^mm/week ,  it 
is  in  the  "static  type"  propagation  region.  The  fracture  morph¬ 
ology  is  predominantly  intergranular  fracture  and  fiber  tear. 

The  entire  da/dN-K  curve  can  be  approximately  described  as  an  "S" 

T1  21 

shaped  curve ,  as  shown  in  Figure  1 .  ' 

In  the  recent  twenty  years,  people  have  conducted 
many  studies  on  the  propogating  behavior  of  fatigue  cracks  in 
the  second  stage  and  have  obtained  some  profound  understanding.  1 
third  stage  of  crack  propagation,  however,  is  not  too  important 
to  the  lifetime  estimation  of  structural  parts  because  of  it3 
rapid  development.  For  this  reason,  the  first  stage  of  crack 
propagating  behavior,  which  has  a  greater  sensitivity  to  the 
structure  and  the  environment,  has  attracted  more  and  more  at¬ 
tention.  Especially  with  the  development  of  fracture  mechanics 
and  application  of  new  testing  techniques,  a  large  number  of 
scientific  papers  have  been  published.  Many  mechanical  models  and 
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microscopic  mechanisms  were  proposed.  However,  mutual  contradic¬ 
tions  often  exist  between  these  data  and  results.  Even  opposite 
conclusions  were  obtained.  For  this  reason,  a  brief  review  of 
the  studies  in  recent  years  on  the  microstructure  of  the  fatigue 
crack  in  its  initial  stage,  its  mechanical  behavior,  and  its  in¬ 
fluencing  factors  is  needed. 


Fif.l 


Schematic  of  da/dN  can*  and  mechanisms  of  fracturecn 
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Figure  1 .  Schematic  of  da/dN  curve  and  mechanisms  of  fracture 
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1.  da/dN  (mm/week) 

2.  lattice  constant/week 

3.  threshold  value  AK. 

4.  final  destruction 

5.  lmm/sec. 

6 .  1mm /hour 

7.  lmm/day 

8 .  lmm/week 

9.  crack  propagation  rate  at  50Hz 

10 .  Segment  A 

11.  non-continuous  mechanism 

12.  heavy  influcencing  factors 

13.  1.  microstructure 

14.  2.  average  stress 

15.  3.  environment 

16.  Segment  B 

17.  continuous  mechanism  (increasing  band) 
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18.  light  influencing  factors 

1?.  1.  micro structure 

20.  2.  thickness 

21.  heavy  influencing  factors 

22.  3.  some  combination  of  environment,  average 
stress,  and  frequency 

23.  Segment  C 

24.  static  mechanism  (cleavage,  intergranular  fiber) 

25.  heavy  influencing  factors 

26.  1.  microstructure 

27.  2.  average  stress 

28 .  3 .  thickness 

29.  little  influencing  factor 

30 .  4 .  environment 

II.  MICROSCOPIC  CHARACTERISTICS  OF  INITIAL  CRACK  PROPAGATION 

[31 

Since  Forsyth  studied  the  fatigue  propagation  of  aluminum 
alloys  in  the  early  sixties  and  divided  it  into  the  I  and  II 
propagating  stages ,  a  great  deal  of  exploratory  studies  of  the 
propagation  behavior  along  the  45°  crystallographic  axis  in  the 
stage  I  of  propagation  have  been  carried  out.  Furthermore,  some 
physical  models  were  presented.  The  details  are  in  Reference  [2] . 

In  the  aspect  of  the  microscopic  morphology  of  fatigue  crack 

r  4  7 

in  the  metal  propagation  period,  Hertzberg  and  Mills  discover¬ 
ed  that  when  the  propagation  rate  was  extremely  low,  the  fracture 
showed  a  crystallographic  facet  or  river  pattern  in  their  obser¬ 
vation  of  the  micro fractography  of  Al ,  Cu,  Ti  and  Ni  alloys,  and 

f  51 

stainless  steel .  Since  then ,  Beevers  summarized  the  work 
done  by  previous  workers.  He  listed  the  crystallographic  direc¬ 
tions  of  the  facets  produced  in  the  initial  propagation  in  over 
ten  alloys.  These  facets  are  usually  close  to  the  densely  packed 
crystallographic  plane  of  the  matrix,  such  as  the  {ill}  and  { 00 if 
planes  of  a  face  centered  cubic  metal,  the  { 1 10 }  and  |00l}  planes 
of  a  body  centered  cubic  metal,  and  f 0001 f  plane  of  a  tightly 
packed  hexagonal  metal .  These  facet  fractures ,  which  emerge 
under  an  extremely  low  propagation  rate,  had  been  believed  to  be 
due  to  the  repeated  slip  in  the  local  region  at  the  tip  of  the 
crack.  This  causes  the  binding  strength  of  atomic  bonds  between 
two  neighboring  slip  surfaces  to  weaken.  Therefore,  the  local 
slip  surface  created  a  cleavage  fracture  under  a  low  tensile 
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stress.  J  However,  the  fact  that  some  of  the  facet  directions 
of  certain  cubic  metals  appeared  on  the  non-slip  surface  {00lj 
was  not  explained. 

However,  when  we  considered  the  crystal  structure,  slip  face 
energy  (SFE) ,  and  modes  of  slip  of  metals  and  alloys,  it  was  dis¬ 


covered  that  the  orientation  of  these  facets  had  a  certain  pat¬ 


tern 
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as  shown  in  Table  1 . 
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//•  x 
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7-  Ni** 

/a.  c 

/S'.  T* 

urn 

/3.  * 

urn 

BCC 

Fc-3XSi 

If  * 

<001> 

?. 

/&  * 

(011>±10* 

HCP 

/O.  Ti6* 

/<$.  X 

fooouts* 

Table  1.  Crystallographic  orientations  of  facets  on  fracture 
surfaces 


1. 

lattice 

2. 

alloy 

3. 

SFE 

4. 

mode  of  slip 

5. 

orientation  of  facet 

6. 

Al  alloy 

7. 

Ni  alloy 

8. 

stainless  steel 

9. 

8  brass 

10. 

Ti  alloy 

From  Table  1  one  can  see  that  it  is  favorable  to  show  a 
cleavage  fracture  along  the  cubic  plane  for  cubic  metals  with 
wavy  slips  and  high  SFE  because  it  was  easy  for  slips  to  occur. 
For  example,  certain  cubic  metals  such  as  Fe-3%Si  could  be  de- 
scribed  by  the  fracture  mechanism  proposed  by  Cottrell1 2 3 4 5 6  ,  i.e. 

~  (III)  + -|- (111)—  a  (001)  (A) 

For  other  cubic  metals  (Ni-based  alloys,  stainless  steel, 

and  B  brass,  etc.),  the  aforementioned  dislocation  is  difficult 

to  occur  due  to  the  low  SFE  and  planar  slip.  Consequently, 

r  7 1 

cleavage  only  emerges  along  the  slip  surface1  .  Recently,  in 

the  analysis  of  the  mechanism  of  the  initial  propagation  of 
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fatigue  crack  by  Lynch1  J ,  it  was  believed  that  the  propagation 
of  crack  along  the  {00l}  plane  in  the  <110>  direction  is  caused 
by  the  alternative  slip  of  the  tip  of  the  crack  in  two  {ill} 
planes,  as  shown  in  Figure  2. 


Figure  2.  Schematic  diagram  of  |l00}  {lio}  propagation  by  .g. 

alternative  slip  in  two  {ill}  planes  at  crack  tip1 

1 .  crack 

2.  slip  plane  1 

3.  slip  plane  2 

4.  (Ill)  slip  plane 

5.  (100)  slip  plane 

6.  (Ill)  slip  plane 


In  practical  engineering  materials,  due  to  the  influence  of 
structure  (the  state  and  distribution  of  the  second  phase  and 
impurities),  environmental  media,  and  stress,  the  microfracto- 
graphic  state  during  the  initial  stage  of  crack  propagation  is 
often  various.  Gerber ich  and  Moody ,  in  their  observation  of 
the  fracture  morphologies  of  structural  steel  and  titanium  alloy, 
have  pointed  out  that  the  microstructure  of  the  fracture  of 
these  two  materials  in  the  initial  propagation  stage  could  be 
divided  into  the  intragranular  and  intergranular  type .  Further¬ 
more,  it  has  nine  patterns:  (a)  inclined  cleavage;  (b)  degener¬ 
ate  cleavage;  (c)  crack  between  a/p  phases;  (d)  a  phase  cleavage 
in  a/p  region;  (e)  intergranular  crack  in  hard  matrix;  (f)  inter¬ 
granular  crack  in  soft  matrix;  (g)  elastic  inte'rgranular  stripe; 

(h)  grain  boundary  void  created  by  point  mass;  and  (i)  linkage 
of  intergranular  crack,  as  shown  in  Figure  3.  The  first  four 
types  are  intragranular  forms  and  the  latter  five  are  intergranu¬ 
lar  forms. 

Freudental^11^  believed  that  because  there  is  always  in¬ 
homogeneity  and  discontinuity  in  the  local  structure  in  a  prac¬ 
tical  structural  part  and  there  are  various  distortions  on  the 
surface  and  in  the  inside,  the  microscopic  mechanism  of  fatigue 
crack  is  primarily  divided  into  the  slip  plane  crack  produced  by 
the  repeated  slip  of  the  soft  metal  and  the  shear  crack  of  high 
strength  alloy  without  slip.  In  the  medium  strength  alloys,  the  /16 
mechanism  is  a  combination  of  both. 
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Figure  3.  Types  of  microfractographs  of  fatigue  crack  propaga¬ 
tion  at  initial  stage 

One  can  see  that  in  practical  engineering  structural  parts 
and  specimens,  the  fracture  structures  are  various.  However, 
they  also  have  a  certain  basic  variation  pattern.  Ouyang  Jie^12^ 
summed  up  eight  types  of  fractures  based  on  the  orientation  of 
cracks  and  fracture  analysis  of  a  large  number  of  specimens  and 
parts  made  of  aluminum  alloys,  structural  steel,  and  high  temper¬ 
ature  super  alloys.  Furthermore,  he  believed  that  these  frac¬ 
tures  could  basically  be  described  by  the  "alternative  shear" 

.  r  2  T 

mechanism  and  the  "alternate  tear-shear"  mechanism.  The  author 

believed  that,  from  a  simple  stress  analysis,  it  belongs  to  the 
dual  directional  tensile  plane  stress  type.  In  this  case,  be¬ 
cause  the  tangential  stresses  on  the  slip  plane  cancel  each  other 
and  the  normal  stress  superimposes  itself,  it  favors  the  occur¬ 
rence  of  a  facet  type  cleavage  fracture. 

III.  MECHANICAL  BEHAVIOR  AND  THRESHOLD  VALUE  IN  INITIAL  CRACK 
PROPAGATION  STAGE 

In  recent  years,  in  order  to  save  energy,  reduce  material 
consumption,  and  improve  product  function,  it  is  unavoidable  to 
have  notches  and  stress  concentrating  zones  present  in  engineer- 
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ing  structural  parts.  Especially  certain  components,  such  as 
rotor  parts  for  an  engine,  gears,  crankshafts,  pressurized  con¬ 
tainers,  etc.,  frequently  sustain  small  amounts  of  low 
frequency  heavy  loads.  In  addition,  the  periodical  inspection 
of  these  parts  in  long  term  operation  is  relatively  difficult. 

For  this  reason,  these  parts  are  required  to  have  an  assurance 
of  safe  crack  propagation  (or  no  propagation).  Consequently,  a 
new  important  branch  has  been  developed  from  "damage  tolerance 
design"  in  product  design,  i.e.,  "safe  crack  design. "  Hence, 

the  further  study  of  mechanical  behavior  and  threshold  values  of 
initial  propagation  of  fatigue  crack  is  seriously  looked  upon. 


1.  Fatigue  Limit  and  Crack  Propagation  Threshold 

Figure  4(a)  is  a  schematic  diagram  of  the  stress-strain  field 

on  a  notched  specimen  and  the  inception  and  propagation  zone  of 

fatigue  crack.  From  the  fatigue  one  can  see  that  because  of  the 

presence  of  a  relatively  large  plastic  and  elastic  pressure 

stress  field  in  front  of  the  tip  of  the  crack,  the  crack  stops 

or  propagates  at  an  extremely  low  rate.  Figure  4(b)  is  the 

schematic  diagram  of  the  relation  between  smooth  and  notched 

r  1 3 1 

fatigue  limits  and  the  fatigue  threshold  stress. 1 

r  1 4 1  ,  ,  . 

Lukas 1  J  pointed  out  that  there  were  three  fatigue  limits 

at  the  present  moment,  i.e. ,  crack  propagation  threshold  (AKfch) , 
fatigue  stress  limit  (<7a)  ,  and  fatigue  plastic  strain  limit 
(eapC).  The  latter  two  fatigue  limits  have  the  following  re¬ 
lation,  i.e.. 


o.-Ke;,, 


(1) 


where  K  and  n  are  material  constants.  The  fatigue  threshold  and 
fatigue  stress  limit  also  have  a  correlation.  Assuming  that 


Katb  is  the  basic  threshold  value,  under  the  cycling  of  a  sym¬ 
metric  stress,  K<*tb  =  AKth^'  For  a  Perf°ratec*  crack,  we  can 
obtain  that 


K,t ,—  1.120.1/  na, 

where  a  is  the  critical  crack  dimension, 
e 


(2) 

When  the  specimen  has 


a  notch 


(3) 

where  oq  is  the  notch  fatigue  limit,  K.  is  the  notch  stress  con-  /17 
n  ^ 

centrating  coefficient,  and  p  is  the  radius  of  the  notch.  The 

aforementioned  formula  can  be  written  approximately  as 


with  increasing  strength,  oG  gradually  increases  and  Katb  de¬ 
creases.  The  second  term  in  the  parenthesis  will  also  decrease. 
Oq  will  be  close  to  the  value  of  (7C/Kt.  When  iKth  is  already 
known,  the  simplified  formula  can  be  applied  in  practical  engi¬ 
neering. 

Smith recently  presented  the  correlation  between  the 
smooth  and  notched  fatigue  limits  and  the  thresholds  in  a  com¬ 
prehensive  analysis  article  regarding  fatigue  limits  and  thres¬ 
hold  values.  According  to  the  usual  fracture  mechanics  formula 
(LEFM) ,  the  threshold  value  (AKth)  is  a  function  of  the  exter¬ 
nally  applied  stress  (to)  and  the  crack  length  (a).  When  a  crack 
begins  to  propagate 

(5) 


where  f  is  the  calibration  coefficient  of  the  geometric  dimension 
when  a  penetrating  crack  is  formed  on  the  surface. 


(6) 


In  order  to  make  the  crack  not  propogate,  it  is  necessary  that 


(7) 

From  the  aforementioned  formula  one  can  see  that  when  the  crack 
length  is  very  small,  it  is  necessary  to  apply  a  larger  stress 
amplitude  to  make  the  crack  propagate.  In  Figure  5,  the  thres¬ 
hold  value  AKth  of  a  long  crack  can  be  expressed  by  a  straight 
line  with  a  slope  of  — *s ,  i.e., 
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240  x 


For  soft  steel  (AKfch  =  9.9  x  106  NM"3/2,  a<tq  = 

10®  NM-*5)  ,  ao  is  approximately  0.2mm.  For  hardened  steel  (AKth= 
3.4  x  10®  NM"3/2,  a<70  =  670  x  10®  NM”*5)  ,  aQ  is  only  0.005mm. 

This  is  why  the  surface  scratches  will  greatly  reduce  the  fatigue 
lifetime.  This  is  because  the  presence  of  scratch  may  reduce  or 
even  eliminate  the  inception  and  initial  propagation  stages  of 
cracks . 


Figure  4. 


(a)  Stress/strain  field  in  a  notched  specimen  and  (b) 
relationship  of  fatigue  limit  of  plain,  notched 
specimens  and  threshold  stress 

1.  <r/e  notch  plastic  region 

2.  crack  tip  region 

3 .  base  material 

4.  notch  elastic  region 

5.  maximum  stress  o>  „ 

max 

6.  safe  band 

7.  plain  fatigue  strength 

8 .  fracture  band 

9.  crack  not  propagating 

10.  stress  concentrating  coefficient  Kt. 
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In  the  meantime,  one  can  also  see  that  when  the  crack 
length  is  less  than  aQ,  some  of  the  formulas  established  in  con¬ 
tinuous  mechanics  (LEFM)  are  no  longer  applicable.  In  this  case, 
the  crack  propagation  pattern  is  different  from  the  one  observed 
with  long  cracks,  as  shown  in  Figure  5(b).  This  phenomenon  was 
also  proven  by  the  experimental  results  recently  published  by 
Romaniv  et  al. When  the  crack  length  exceeded  its  critical 
value  (ag°)  ,  aK^  was  not  related  to  the  crack  length.  They  also  /18 

gave  the  correlating  formula  between  the  short  crack  threshold 
sh 

value  (aK  h)  and  the  long  crack  threshold  value  (aK  ) 


where  r  is  a  material  constant. 


Figure  5.  (a)  Transition  from  fatigue  limit  of  a  plain  specimen 

to  K..  of  a  long  crack  and  (b)  its  effect  on  the 
measurement  of  threshold  valued 15 1 

1.  ^  oQ  constant 

2.  Propagation  Behavior  of  Short  Cracks 

Kitagawa  and  co-author  have  studied  the  threshold  stress 
of  structural  steel  rod  specimens  of  various  strength  levels  with 
microcracks  on  the  surface,  and  obtained  the  following  correlation 


A  o,k 


A  K,h 

f  (  a  )  V  xa 


(10) 


where  f(a)  is  a  function  related  to  the  crack  size  and  specimen 
geometry . 

El  Haddad  et  al^^'^1  also  further  verified  this  phenome¬ 
non  using  a  notched  thin  plate  CSAG40-11  steel  specimen.  In  order 
to  make  the  method  of  LEFM  fracture  mechanics  effective,  they 
introduced  a  constant  lo  for  correction,  i.e., 

Atf-A^a  (a  +/,)F  (ID 

where  F  is  the  coefficient  related  to  the  geometrical  size  of  the 
specimen.  For  a  small  penetrating  crack,  F  is  approximately  equal 
to  1.  When  a  is  approximately  zero,  lQ»a,  tjj=  o L,  and  AK  = 

AKth*  We  °btaine<^  a  correlation  similar  to  equation  (8)  . 


, _ 1  /  A/C,.  V 

*  *  V  Ao,  J 


(12) 


When  considering  a  crack  near  the  surface  using  1.12  as  the 
coefficient,  equation  (8)  is  basically  in  agreement  with  (12). 

By  substituting  equation  (11)  into  the  above  formula,  we  can 
obtain  the  general  formula  for  the  threshold  stress  (Acr_h)  at 
any  crack  length. 


Ao,.- 


A  A,. 


(13) 


Figure  6  shows  the  curve  correlating  the  fatigue  threshold 
stress  and  the  crack  length  of  G40-11  steel.  One  can  see 

that  when  tQ  =  0.24mm,  the  estimated  value  agrees  with  the  ex¬ 
perimental  results  well.  With  regard  to  the  physical  signifi¬ 
cance  of  1Q,  it  is  not  clear  at  the  present  moment.  It  may  be 

related  to  the  stress-strain  field  at  the  tip  of  the  crack. 

r  i  g  I 

Topper  and  El  Hadded1  J  obtained  a  similar  formula  for  the  /19 
propagation  of  a  short  crack  through  an  analysis  of  the  initial 
propagation  of  fatigue  cracks  on  notched  specimens . 


Aor,. 


a/c,. 

V  *  ( /  +/0) 


(14) 


where  K1  is  the  elastic  stress  concentrating  coefficient,  and  Z 
is  the  crack  length.  The  above  formula  can  be  written  as 


For  the  center  or  edge  notch  with  a  depth  c 


When  I  approaches  zero,  K’  is  close  to  1.12  K 


(16 


to' 


From  equations  (15)  and  (16)  one  can  get  that 


rx  /  /. 


V  /  +  c\  i+i. 


(r 


2. 

MQ 


Figure  6.  Fatigue  threshold  stress  versus  crack  length  for  a 
G40-11  steel  l17l 

1.  a  (mm) 

2 .  A  <r  “constant 

3.  I  =0 .24mm 

4.  experimental  value  of  G40-11  steel 

5 .  calculated  value 

6.  Aa../2  (106N.M“1 2 3 4 5 6) 


Figure  7.  Comparison  of  calculated  curves  with  the  experimental 

data  118] 

1 .  £  (mm) 

2.  V-shaped  notch  c=5m 

3 .  calculated  value 

4.  1=0 .26mm  .e. 

5.  experimental  results1  J 

6 .  fracture 

7.  0.25mm 

8.  1.27mm 

9 .  7 . 6mm 

10.  crack  not  propagating 

11.  ififf/ 2  (10®  N.M-2) 


Act, 


Fact. 

1  +  V  c/l. 


(18) 


Figure  7  shows  the  comparison  of  the  experimental  values  with 

calculated  values  obtained  using  equations  (15)  and  (18)  for 

specimens  with  various  notch  diameters  (p)  (V  shaped  notch,  depth 

c=5mm,  £  =0 . 29mm) . 
o 

For  a  "sharp  notch,"  the  crack  inception  stress  is 


Ao,-Aa,//f, 


(19) 


This  stress  increases  as  £  extends  until  reaching  its  maximum  /20 
value . 
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For  a  "dull  notch",  this  stress  decreases  with  increasing 
!(or  equal  to  the  initial  stress),  i.e. 

Kt^l  +  2v/c/P  (20) 

The  boundary  line  of  the  two  types  of  notch  is  P  =4i  .  From 

cr  ® 

equation  20  we  can  obtain 

I  +  2  /c/P,,  (21) 

Then 

l  yit  (22) 

Figure  8  shows  the  initial  short  crack  propagation  behavior 

r  181 

of  notches  of  various  geometrical  shapes . 

Hudak  Jr.^^l  had  given  the  following  preliminary  conclusions 
through  an  analysis  of  short  crack  propagation  behavior: 

(1)  Using  the  general  fracture  mechanical  analysis  method 
(LEFM) ,  the  propagation  behavior  of  a  short  crack  is  different 
from  that  of  the  long  crack.  The  primary  reason  is  that  the 
assumption  in  the  conventional  continuous  mechanics  is  already  no 
longer  applicable. 

(2)  A  short  crack  might  propagate  below  the  crack  size  es¬ 
timated  and  stress  level  applied  for  a  long  crack. 

(3)  The  growth  rate  of  a  short  crack  is  higher  than  that  of 
the  growth  rate  of  a  long  crack  predicted  by  the  data  da/dN-AK. 

(4)  The  actual  propagation  behavior  of  a  short  crack  should 
take  into  account  of  a  constant  length  t  on  top  of  the  actual 
crack  size.  At  the  present  moment,  the  physical  meaning  of  lQ  is 
not  understood.  However,  it  can  be  used  to  adequately  measure 
the  threshold  value  of  a  short  crack  and  its  propagation  behavior. 


da/JN  (**/«) 


AJC(***-*-*'*> 


30  40  60  80  100 


Figure  8.  Propagation  behaviors  of  short  crack  at  notches 
different  geometry  £ 1 8 ’ 

1.  AK  (101 2 3 4 5 6 7 8N.m"3/2) 

2 .  center  crack 

3.  circular  hole  P/b=0.134 

4.  circular  hole  P/b=0.227 

5.  elliptical  hole  p=1.19 

6.  A<r=269  x  106  N.M-2,  R=-l 

7.  G40-11  steel 

8 .  da/dN (mm/week) 


ttnattSL’S- 


Figure  9.  Elastic  and  plastic  stress-strain  field  at  the  root 
of  a  notch  t 2 1 1 

1 .  notch  plastic  region 

2.  notch  elastic  region 
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Miller  and  his  co-authors  '  had  questioned  the  afore¬ 
mentioned  analysis.  They  believed  that  the  physical  meaning  of 
introducing  the  constant  lQ  was  not  clear .  The  various  propaga¬ 
tion  behaviors  of  some  short  cracks  below  Act  could  not  be  pre- 

e 

dieted  using  equation  (11).  In  addition,  tQ  could  not  have  been 
a  constant  because  it  was  related  to  the  level  of  applied  stress. 
For  this  reason,  they  proposed  a  model  using  a  local  plastic 
shear  model  to  describes  the  inception  and  initial  propagation 
of  a  crack,  as  shown  in  Figure  9.  In  the  figure,  AS  is  the 
plastic  tangential  displacement  of  the  notch.  It  is  a  function 
of  stress  level,  cycle,  yield  strength,  and  notch  size  and  shape. 

When  a  crack  enters  the  plastic  region,  the  shear  displace¬ 
ment  |0C}  of  the  crack  tip  can  be  considered  to  originate  from 
the  plastic  tangential  displacement  {*0p|  caused  by  the  notch 
plastic  region.  This  tangential  displacement  decreases  as  the 
crack  penetrates  into  the  plajstic  region.  In  addition,  the 
plastic  distortion  produced  by  the  strangeness  of  the  crack  tip  /21 
(0e)  must  be  added  (which  can  be  obtained  according  to  the  LEFM 


36 


analysis) .  Therefore 


<j><“4>,+4>. 


(23) 


They  believe  that  the  critical  value  of  0t  for  the  inception  and 
propagation  of  a  crack  could  be 


(24) 


where  is  the  cyclic  yield  strength.  When  a  crack  begins  to 

grow  from  a  notch,  the  rate  of  decrease  of  0p  might  be  faster 
than  0e.  Consequently,  the  situation  that  0t<^th  will  appear,  as 
shown  in  the  1 0  |e2  curve  in  Figure  10.  Thus,  a  non-propagating 
crack  is  generated. 

In  summary,  in  regard  to  the  propagation  and  mechanism  of 
a  short  fatigue  crack  in  the  initial  stage,  it  is  still  in  the 
developmental  stage.  It  is  awaiting  further  experimental  work 
and  more  profound  theoretical  analysis. 


Figure  10. 


Schematic  of  fatigue  crack  propagation  behavior  at  a 
notch  under  various  stress  levels [21] 
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IV.  SOME  FACTORS  AFFECTING  INITIAL  PROPAGATION  OF  CRACKS 


1.  Effect  of  Stress  Ratio 

Among  the  factors  affecting  the  initial  crack  propagation 
behavior  and  threshold  value,  it  is  generally  believed  that  the 

r  2  2  T 

stress  ratio  R  is  one  of  the  important  factors .  The  increase 

in  R  will  cause  the  decrease  of  the  initial  propagation  rate  and 
threshold  value . 

With  regard  to  the  effect  of  the  ratio  R  or  the  average 

stress  on  the  AK  .  value  in  a  relatively  low  propagation  rate 

cn  r  ?  3 1 

range,  it  is  usually  expressed  by  the  following  formula. L  J 


aa:,*=a.k:,4|(  i  -ry 


(25) 


where  r  is  a  material  constant 
when  R=0. 


AK 


thQ 


is  the  threshold  value 


It  is  obvious  that  the  above  formula  is  only  applic- 

(241 

able  to  the  range  where  R  is  positive.  Kaisand  and  Mowbray1 


had  given  a  generalized  formula  which  agrees  with  the  experimen¬ 
tal  data  well , 


&K,k=AK,ktf(R) 

(26) 

when 

R  >0 

(A) 

when 

R  <0 

"“-(or 

(B) 

That  means  when  R  is  negative,  they  assumed  that  it  was  in  a 

loaded  cycle.  Only  the  R  values  in  the  1/3  compression  portion 

are  effective  to  the  propogation  of  cracks. 

There  have  been  some  different  explanations  regarding  the 

effect  of  the  stress  ratio  R  on  the  da/dN  of  the  initial  crack 

r  251 

and  the  threshold  value.  Klesnil  and  Lukas  believed  that 
the  load  experience  is  an  important  factor  for  structural  steel. 
They  also  pointed  out  that  the  effect  of  R  on  the  threshold  value 
is  caused  by  the  compressive  stress  in  the  plastic  region  on  the 
tip  of  the  crack.  Schmidt  and  Paris1  ,  however,  adopted  the 
crack  closure  effect  to  explain  the  effect  of  R  on  the  threshold 


value . 


For  a  low  stress  value,  K  .  «  K  _ 

mm  cl 


K__„  =  K  _  +  AKa  =  constant 
max  cl  0 


AK.  =  K  -K  .  =  K  (1  -  R) 

th  max  mm  max ' 


For  a  high  stress  ratio,  K  .  %  K 


'min*  cl 


AKfch  =  AKq  =  constant 


fC  n  A/C, 

“  l-R  ~l-R 


The  above  relations  have  been  verified  in  2024-T3  aluminum 
alloy,  A533  steel,  and  T-l  titanium  alloy.  Experimental  results 
pointed  out  that  the  load  frequency  has  a  certain  effect  on  the 
aluminum  alloy  crack  closure  stress  strength  factor  KCI  and  AKth. 
At  a  low  stress  ratio,  the  effect  of  frequency  on  K__„  is  not 

r  2  7 1  max 

too  large.  J 


Radon 


had  discovered  that  the  effect  of  R  was  related 


to  the  thickness  of  the  specimen  in  his  study  of  BS-4360  low  alloy 
steel  at  near-threshold  range.  For  a  relatively  thin  specimen 
(B  =*  12mm)  ,  the  effect  was  more  apparent.  In  a  50mm  thick  speci¬ 
men,  however,  the  effect  of  R  was  not  too  large.  This  phenomenon 
might  be  related  to  the  stress  state  of  the  specimen  and  the  clo¬ 
sure  effect. 

In  addition,  the  prestrain  also  has  different  effects  on  the 

threshold  values  of  carbon  steel  at  various  R  ratios  (see  Figure 

(291 

11).  Blacktop  et  al 1  ‘  believed  that  at  a  low  stress  ratio  the 

effect  of  cold  work  might  be  related  to  the  size  of  the  plastic 
zone  of  the  crack  tip.  Simultaneously,  they  also  pointed  out 
that  according  to  the  following  formula 


A/C  2*P,u, 
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<  V. 


V  --  -V.  A 


the  prestrain  will  cause  the  ef  in  the  equation  to  decrease. 
Consequently,  the  threshold  value  will  decrease. 


Figure  11.  Effect  of  prestrain  on  AK  .  at  different  stress r2g, 
ratios,  R,  in  a  plain  careon  steel  (0.1C-1.5Mn) 1  J 

1.  AKth  (106N.M-3/2) 

2.  strain  rate  {%) 

2.  Effect  of  Microstructure 

With  regard  to  the  effect  of  grain  size  on  the  initial 
crack  propagation  behavior  and  threshold  value,  a  great  deal  of 
experimental  data  and  many  empirical  formulas  have  been  publish¬ 
ed  to  date.  For  example,  Kitagawa  et  al*-31^  had  discovered  that 
and  <7y  had  the  following  relation  in  testing  a  low  strength 
alloy  steel 

A  K,»~AE—tiof  (32) 

where  A  and  B  are  experimental  constants.  Beevers^  proposed 
the  following  empirical  formula  in  his  study  of  the  effect  of 
grain  size  (d)  on  the  threshold  value  of  a  manganese  steel 


40 


A/r,*-3.8  +  l.4X10,rf,/1 


0,-139  +  0. 347c<- 


Taira  et  al 1  1  pointed  out  that  AKfch  increased  with  in¬ 

creasing  grain  size  (d*5)  in  their  study  of  the  effect  of  grain 
size  on  the  threshold  value  of  a  low  carbon  steel .  Furthermore , 
a  computation  formula  was  given  for  the  effective  threshold 
(K  £fth)  and  grain  size  d  . 


A#,//,*— 0.19+1.96*  10*  V  d 


(35)  /23 


where  AK  *  K  -  K  .  They  discovered  that  the  crack  tip 
eff  max  op 

slip  belt  region  e  increases  with  the  grain  size  (d) ,  i.e., 

s 

e  =  0.51d.  The  da/dN-AK  curve  of  various  grain  sizes  approaches 
s 

a  straight  line  if  it  is  expressed  by  AK eff,  as  shown  in  Figure 

12. 

Gerberich  and  Moody also  pointed  out  that  the  slope  of 
the  AK^/E-d1*  curve  of  different  materials  showed  a  positive  and 
a  negative  type  through  an  analysis  of  the  experimental  data  of 
structural  steels  and  titanium  alloys  of  different  strength 
levels,  as  shown  in  Figure  13. 
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Figure  12.  The  da/dN-AK  curve  for  low  carbon  steel  with  various 
grain  size!32! 

1 .  material 

2 .  (micron) 

3.  crack  propagation  rate  da/dN  (m/week)  _3/2 

4.  stress  strength  range  AK,  AK  (10°  N.M  /z) 
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Figure  13.  AK^/E  versus  d*5  curves  of  various  materials 

1.  AKth/E (xl05m^) 

2.  low  strength  steel,  r  =  0.79 

3.  hjgh  strength  steel,  r  =  -0.53 

4.  d*5  (micron**) 

From  the  results  mentioned  above,  one  finds  out  that  it  is 
not  enough  to  only  consider  the  grain  size.  It  is  also  necessary 
to  analyze  the  internal  microstructure  and  morphology.  Gu  Hai- 
cheng  and  Knott gave  the  following  correlation 

(36) 

They  believed  that  the  threshold  value  is  not  only  related  to 
the  tensile  strength  (Oy)  of  the  material,  but  also  correlated 
to  the  "spacing"  parameter  (p*)  of  the  structure  (i.e.,  Neuber 
material  constant) .  The  increase  in  the  tensile  strength  fre¬ 
quently  brings  about  a  decrease  in  the  structural  spacing  par¬ 
ameter.  Conseugently ,  certain  types  of  steel  may  show  a  peak 
value  of  AK..  at  different  temper  temperatures. 


As  for  the  multi-phase  metallographic  aspect,  Masounave  et 
[34] 

al  studied  the  fatigue  crack  propagation  behaviors  of  various 

ferrites  (a)  and  pearlites  (p)  in  carbon  steel.  They  pointed  out 
that  the  following  relation  between  and  the  microstructure 

AAT,4t-/:<t.+  (i  -/:H,  (37) 

where  fa  is  the  percentage  of  ferrite  by  volume.  0a  and  0p 
represent  the  action  function  of  AK^  with  respect  to  a  and  p. 

0a  =  AKq  +  Kfd^  (38) 

0p  -  AK  =  constant  (39) 

[351 

The  author  and  co-workers  discovered  that  the  functional 
relation  between  AKfch  and  the  microstructure  can  be  approximately 
expressed  by  the  following  empirical  formula  in  the  determination 
of  a  threshold  value  of  a  high  strength  steel  after  different 
temperature  treatments: 

Aif,»«1.95/ii+7.53/#+14.l/*  (40) 

where  fM,  fn,  and  f.  are  the  volumetric  percentages  of  tempered 
martensite,  banite  and  residual  austenite,  respectively.  From 
the  above  formula  one  can  see  that  the  contribution  of  different 
microstructures  to  AKfch  is  not  the  same,  i.e.,  M:B:A  is  approxi¬ 
mately  equal  to  1:4:7.  This  estimation  is  approximately  equiva¬ 
lent  to  the  contribution  of  the  10%  residual  austenite  in  4340 
steel  to  the  plastic  distortion  energy  studied  by  Schwalbe . 

In  addition,  through  the  determination  of  the  lattice  distortion, 
dislocation  density,  residual  stress,  and  microscopic  hardness 
of  the  surface  layer  structure  also  significantly  affected  the 
fatigue  crack  propagation  and  the  threshold  value. (35,37] 

Brown  and  Smith discovered  that  the  corresponding  re¬ 
lation  between  AKfch  and  s/G&  in  a  equi-axed  structure  alloys  with 
residual  8  phase  was  not  as  apparent  as  that  of  a  low  strength 
steel  in  studying  the  effect  of  microstructure  and  crystal 
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orientation  on  the  threshold  of  Ti-6A1-4V.  However,  the  specimen 
T-L  perpendicular  to  the  (0001)  surface  and  the  specimen  L-T  par¬ 
allel  to  the  (0001)  surface  have  different  initial  crack  propaga- 

6  —3 /2 

tion  behaviors.  For  the  latter,  when  AK*3.5  x  10  Nm  ,  a 
crack  propagation  threshold  value  has  still  not  been  found. 

Looking  at  the  experimental  results  mentioned  above ,  the 
effect  of  the  microstructure  on  the  initial  crack  propagation  and 
the  threshold  value  is  very  complicated.  Some  explanations  have 
been  given.  In  general,  it  is  believed  to  be  closely  related  to 
the  reverse  plastic  zone  size  (ry) , 


The  relative  dimension  between  ry  and  the  microstructure  is  an 
important  factor  affecting  the  initial  crack  propagation  and 
threshold. 

[391 

Irving  and  Beevers  believed  that  the  growth  of  a  crack 
became  structure  sensitive  from  structure  insensitive  as  the 
size  of  the  plastic  zone  reached  the  microstructure  dimension 
when  they  studied  the  effect  of  the  microstructure  of  the  Ti- 
6A1-4V  titanium  alloy  on  fatigue  crack  propagation.  Kao  and 
Byrne also  pointed  out  that  the  cementite  spacing  in  pearlite 
steel  is  the  smallest  microstructure  unit.  It  has  a  primary  in¬ 
duction  effect  on  the  fatigue  crack  propagation  of  steel.  The 
effect  of  the  original  austenite  grain  size  and  its  grain  bound¬ 
ary  on  the  crack  propagation  in  the  near  threshold  range  belongs 
to  a  secondary  type .  Furthermore ,  the  spacing  determines  the 
free  mean  distance  of  the  crack  leading  edge  propagation  dis¬ 
placement.  A  decrease  in  the  spacing  limits  the  displacement  of 
the  leading  edge  of  a  crack.  Consequently,  its  propagation  rate 
is  reduced.  This  effect  is  even  more  apparent  at  the  near¬ 
threshold  range. 

3.  Effect  of  Environmental  Media 

Among  all  the  influencing  factors,  the  effect  of  environ¬ 
mental  media  is  a  subject  which  attracts  a  lot  of  attention 
presently  and  is  very  complex.  In  the  past,  because  of  the  lack 


of  experimental  results  in  the  initial  crack  propagation  stage , 
some  mechanisms  in  the  middle  and  late  stages  of  crack  propaga¬ 
tion  (such  as  metal  dissolution  and  hydrogen  embrittlement)  were 
frequently  used  to  explain  it.  Recently,  Ritchie ^  summarized 
the  following  points  according  to  the  experimental  results  ob¬ 
tained  in  the  aspect  of  the  effect  of  environmental  media  on  the 
behaviors  in  the  initial  crack  propagation  period  and  in  the 
threshold  range: 

(1)  the  effect  primarily  occurs  in  a  low  stress  ratio 
condition 

(2)  as  compared  with  the  results  obtained  in  laboratory  air, 
the  dry  inert  atmosphere  accelerates  da/dN  and  it  decelerates  in 
water . 

(3)  for  a  high  strength  steel,  the  propagation  in  dry  hydro¬ 
gen  is  lower  than  that  in  laboratory  air. 

Figure  14  shows  the  propagation  behaviors  of  2k  Cr-lMo  bani- 

tic  steel  in  moist  air  and  dry  hydrogen  environment  at  various 

stress  ratios  and  frequencies.  One  can  see  that  the  effect  of 

hydrogen  on  crack  propagation  appeared  in  two  regions,  i.e.,  at 

-5 

a  medium  rate  (da/dN~10  .mm/week)  and  around  the  threshold  value 
(da/dN<10~^mm/week) 

The  effect  of  the  environment  on  the  crack  propagation  of 
300M  high  strength  steel  is  slightly  different,  as  shown  in 
Figure  15.  When  the  rate  is  higher  than  10  ^mm/week,  hydrogen 

g 

causes  an  acceleration  (f=5Hz) .  When  it  is  lower  than  10  mm/week 

T431 

dry  hydrogen  reduces  the  mitral  crack  propagation  rate  by  16% i  J 

Such  a  decrease  is  also  shown  in  other  high  strength  rotor 
r  44 1 

steels.  In  the  salt  mist  test  of  high  strength  steel,  such 

a  decreasing  effect  was  also  observed.  The  effect  of  dry  hydro¬ 
gen  on  the  crack  propagation  behavior  of  structural  steel  is 
schematically  shown  in  Figure  16. 

The  mechanism  concerning  the  effect  of  the  environmental 
media  on  the  initial  crack  propagation  could  be  approximately 
divided  into  two  main  groups,  i.e.,  the  hydrogen  embrittlement 
model  and  the  crack  closure  model.  Although  the  former  success¬ 
fully  described  the  behavior  of  a  fatigue  crack  at  a  medium  or 


relatively  high  propagation  rate ,  yet  according  to  the  experi¬ 
mental  results  described  above  there 


Figure  14.  Fatigue  crack  propagation  behavior  in  moist  air  and 
dry  hydrogen  for  a  banitic  SA542-3  steel 

1.  crack  propagation  rate  da/dN  (mm/week) 

2.  2\  Cr-lMo  steel 

3.  SA542-3  frequency  50Hz,  room  temperature 

4 .  moist  air 

5 .  dry  air 

6.  (0.5  Hertz) 

7 .  air 

8 .  hydrogen 

9 .  hydrogen 

10 .  air 

11.  crystal  lattice  constant/week 

12.  threshold  value  AKfl  ,  _ 

13.  stress  strength  range  AK(10°  N.M~  /  ) 
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Figure  16 . 


Schematic  diagram  of  the  effect  of  dry  gaseous  hydro 
gen  on  .fatigue  crack  propagation  in  low  strength 
steels  L J 


1.  da/dN  (mm/week) 

2 .  increasing  R 

3 .  hydrogen 

4 .  increasing  R 

5 .  decreasing  r 

6 .  air  (not  related  to  r  and  R) 

7 .  air  (not  related  to  r 

8 .  AKn  (hydrogen) 

9.  aKq  (air) 

is  obviously  another  mechanism  present.  For  this  reason,  the 
emphasis  is  placed  on  the  crack  closure  model. 

In  the  recent  two  years,  Ritchie  and  co-worker ^ pre¬ 
sented  two  models  for  the  crack  closure  effect,  i.e.,  the  "oxide 
induced"  and  "roughness  induced"  closure  models,  based  on  the 
great  deal  of  experimental  results  obtained  for  the  effect  of 
environmental  media  on  alloy  steels  and  high  strength  steels. 

The  schematic  diagrams  for  the  two  models  are  shown  in  Figure  17 
The  former  is  the  foreclosure  caused  by  the  accumulation  of 
oxides,  while  the  latter  is  caused  by  the  roughness  of  the  frac¬ 
ture  together  with  the  second  type  shear  displacement.  Appar¬ 
ently,  the  closure  caused  by  plastic  distortion  crushes  the 


oxide  film  and  creates  some  fresh  surface.  Consequently,  the 
oxide  film  is  thickened.  They  found  that  the  oxide  film  of 

_3 

4?s  Cr-2Mo  steel  near  the  threshold  could  reach  0.2mm  ,  which  is 

approximately  twenty  times  higher  than  the  thickness  in  the 
medium  and  high  propagation  rate  regions.  The  thickness  of  the 

O 

corrosion  product  in  dry  hydrogen,  however,  is  only  5Q'~-,150A. 

At  a  high  stress  ratio.  Type  II  displacement  will  cause  "corro¬ 
sion-oxidation"  to  significantly  reduce  the  thickness  of  the  oxide 
layer.  Figure  18  shows  the  variation  curves  of  the  oxide  layer 
thicknesses  and  da/dN  of  the  two  types  of  steel  mentioned  above 
versus  the  crack  length  as  measured  by  an  Auger  Spectrometer. 
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Figure  17.  Schematic  of  the  mechanism  for  near-threshold 
fatigue  crack  propagation^ 

1.  no  closure  effect 

2.  oxide  induction 

3.  roughness  induction 


Figure  18.  Auger  measurements  of  crack  flank  oxide  thickness  in 
2  Cr-lMo  and  300-M  steels  as  a  function  of  crack 
length  and  da/dN^2^ 

1.  da/dN  (mm/ week) 

2.  oxide  layer  thickness  (micron) 

3.  air  R=0.05 

4 .  hydrogen  R=0 . 0  5 

5.  air  R=0.75 

6.  air,  hydrogen,  R=0.05 

7.  threshold  AKq 

8 .  crack  length u (mm) 

f  451 

Suresh  et  al  also  derived  a  calculation  formula  for  the 

crack  closure  stress  strength  factor  Kc^ 

K  I  JE 

-o')  (42) 

where  d  is  the  maximum  thickness  of  the  oxide  layer,  and  2j£  is 

the  distance  from  the  maximum  thickness  point  to  the  tip  of  the 
2 

crack.  E/(l-u  )  is  the  effective  elastic  modulus  under  a  plane 
strain  condition.  The  calculated  values  obtained  using  equation 
(42)  are  in  approximate  agreement  with  the  experimental  values. 


Figure  19  shows  the  1:1  corresponding  relation  between  the 
maximum  oxide  layer  thickness  and  the  crack  tip  opening  displace¬ 
ment  ( CTOD )  plotted  using  some  experimental  data  of  2\  Cr-lMo 
steel  in  the  atmosphere,  hydrogen,  helium,  and  water.  A  large 
microstructure  or  a  high  contaminant  content  may  produce  a  rough 
or  facet  fracture  surface  to  promote  the  occurrence  of  the 
"roughness  induced"  closure.  This  result  agrees  with  the  obser¬ 
vation^6-^  that  a  large  '’mismatch"  existed  in  th-i  fracture  sur¬ 
face  of  a  large  a  titanium  structure  in  the  crack  propagation 
near  the  threshold  as  reported  by  early  workers . 

The  above  was  a  discussion  on  some  major  factors  concerning 
the  initial  crack  propagation  behavior.  •  In  reality,  other  factors 
such  as  load  history,  frequency,  amplitude,  load  type,  and  test-  /27 
ing  method,  as  well  as  the  combined  effect  of  various  factors, 
must  be  seriously  considered.  Due  to  the  limitation  in  page 
number,  they  cannot  be  individually  listed.  Some  of  the  recently 
published  overall  review  literature  such  as  [1,  2,  47,  48]  should 
be  referred  to. 
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Figure  19.  Showing  the  1:1  correspondence  of  maximum  oxide 

thickness  with  crack  tip  opening  displacement  (CTOD) 
for  2JsCr-lMo  steel  in  air,  hydrogen,  helium  and 
water lU 

1.  crack  tip  opening  displacement  (micron) 

2 .  2h  Cr-lMo  steel 

3.  maximum  oxide  layer  thickness  d (micron) 

V.  CONCLUSIONS 

In  this  paper,  a  brief  review  was  carried  out  regarding  the 
results  of  recent  studies  of  metal  fatigue  crack  propagation  be¬ 
haviors  and  fatigue  threshold  values.  Some  analysis  and  discus¬ 
sion  were  made  on  the  microscopic  characteristics ,  mechanical 
behavior,  and  its  major  affecting  factors  of  the  initial  crack 
propagation  process .  Because  crack  propagation  in  the  near 
threshold  range  is  structure  sensitive,  also  because  of  the  sig¬ 
nificant  effect  of  the  environment  (media,  temperature) ,  the  load 
condition  and  the  testing  method,  it  is  necessary  to  seriously 
distinguish  and  analyze  the  experimental  results  and  data  pub¬ 
lished.  For  example,  the  humidity  in  the  atmosphere  has  a  strong 
effect  on  certain  materials  (such  as  structural  steel) .  In  the 
determination  of  the  initial  crack  propagation  rate  and  the 
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threshold  and  its  application  in  engineering,  we  must  seriously 
consider  some  of  the  influencing  factors  and  the  combined  effect. 
In  the  area  of  experimental  plan  design,  it  is  necessary  to  suf¬ 
ficiently  combine  the  operating  conditions  of  the  products  such 
as  the  environment,  temperature,  stress  ratio,  etc.  to  find  the 
characteristics  and  patterns  under  specific  conditions.  In  the 
meantime,  it  is  also  necessary  to  combine  the  microstructure  with 
the  macrostructure,  and  to  combine  the  materials,  techniques,  to 
the  structure,  in  order  to  find  the  general  pattern,  to  develop 
the  application  of  fracture  mechanics ,  and  to  explore  ways  to  im¬ 
prove  the  fatigue  characteristics  of  materials  and  structural 
parts.  It  can  then  be  made  to  more  effectively  guide  product 
design,  scientific  studies,  and  production  to  better  serve  the 
economic  construction  for  the  people . 

During  the  organization  process  of  this  report,  assistance 
was  received  from  comrades  Liu  Caimu,  Zhang  Shijie,  Ouyang  Jie, 
and  Yuan  Gaoming,  as  well  as  the  people  in  the  information  office 
I  want  to  express  my  gratitude  here. 
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Finite  Element  Models  Based  on  a  Generalized  Hamilton's  Principle 
for  Pfluger's  Rod  and  Its  Related  Plate* 
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ABSTRACT 

✓ 

A  generalized  Hamilton's  principle  for  analysis  of  Ijluger's 
rod  and  its  relative  plate  is  provided  in  engineering  form  by  a 
quite  simple  approach,  so  that  the  problem  of  a  nonconservative 
system  is  transformed  into  that  of  a  generalized  conservative 
system.  According  to  this  variational  principle,  an  attempt  to 
undertake  the  discrete  tests  of  finite  elements  has  been  made 
and  the  features  of  the  finite  element  model  are  discussed  in 
brief.  A  set  of  symmetric  matrices  for  an  element  model  are  given 
in  explicit  form.  After  it  is  confirmed  that  instability  of  Pfluger's 
rod  is  caused  by  divergence,  not  by  flutter,  it  is  pointed  out  that 
the  critical  load  can  be  determined  by  applying  the  conventional 
iterative  technique  for  Rayleigh  quotient.  Some  numerical 
examples  have  verified  this  model. 


Finally,  as  an  interesting  extension,  a  finite  strip  model 
for  analysis  of  simply  supported  rectangular  plate  subjected  to  follower 
forces  has  been  derived  by  using  the  matrices  mentioned  above  and 
tested  preliminarily. 
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I .  INTRODUCTION 

As  is  well  known,  there  may  exist  in  non-conservative  forces 
so-called  follower  forces  of  the  non-conservative  component.  The 
work  done  by  these  forces  under  non-inf initesimal  displacements  is 
not  determined  by  the  initial  and  final  states  but  depends  on  the 
path  taken.  Hence,  the  problem  of  stability  for  structures  under  the 
action  of  follower  forces,  such  as  the  Pfluger's  rod,  cannot  be  solved 
by  obtaining  an  approximate  solution  in  the  form  of  an  expansion  start¬ 
ing  out  from  the  usual  Hamilton's  principle.  This  has  been  an  establish 
point  of  discussion  for  a  long  time.  Therefore,  if  one  still  wants  to 
handle  these  non-conservative  problems  by  means  of  a  variational 
principle,  one  is  left  with  two  choices.  One  is  to  adopt  a  variational 
principle  that  simultaneously  takes  into  account  the  adjoint  system 
^ .  In  this  case,  the  scale  of  the  problem  to  be  solved  is  enlarged 
because  of  the  introduction  of  the  variables  of  the  adjoint  system. 

The  other  is  to  use  an  incomplete  variational  principle  (or  variational 

f  2 1 

statement)  in  the  expansion  of  an  approximate  solution  .  Although  the 
latter  method  is  a  more  direct  approach,  its  lack  of  symmetry  makes 
numerical  calculations  more  elaborate. 

Is  it  possible  to  establish  a  general  variational  function 
from  a  generalized  Hamiltonian  principle  without  introducing  the 
adjoint  system  and  yet  retaining  symmetry?  This  possibility  exists 
for  problems  with  finitude.  In  Reference  [3],  a  general  discussion 
has  been  given  of  this  method.  However,  its  generality  is  only 
in  form.  In  this  paper,  the  general  variational  function  has 
been  cast  into  a  practical  engineering  form  by  a  rather  simple 
approach  for  Pfluger's  rod  and  a  related  problem.  In  essence,  this 
is  equivalent  to  what  is  given  in  Reference  [3] .  Based  on  this 
variational  principle,  we  have  undertaken  a  first  attempt  at 
finite  element  discretization  tests,  discussed  their  characteristics, 
and  provided  an  element  model  in  explicit  form  whose  matrices  are 
all  symmetrical.  Computed  example  demonstrates  that  this  theory 
also  shows  effectiveness  in  the  numerical  aspect. 


For  a  simply  supported  rectangular  plate  subjected  to  tan¬ 
gentially  distributed  follower  forces,  one  can  use  the  related 
matrices  of  the  Pfluger  rod  element  model  to  obtain  the  corresponding 
matrices  of  the  finite  strip  model  for  the  plate. 

II.  GENERALIZED  HAMILTONIAN  PRINCIPLE  AND  GENERALIZED  LAW  OF 
CONSERVATION  OF  ENERGY 


It  is  easy  to  write  down  the  differential  equation  of  motion 

and  boundary  condition  for  Pfluger' s  rod  as  follows: 

D.  E.i 


B.C.i 


*(*,  i  )  +  -£-»«*(*,  t)+-~-(L~  x)w„(x,  0=0 

ui(  0,  t)—ui(L,  i  )  —  u>„(  0  ,  t )  =*u>n(Z.,  0=0 


(1) 

(2) 


In  the  above  equations,  a  (x,t)  denotes  deflection,  a  is 
stiffness  to  bending,  g  is  the  follower  force  of  the  distribution  in 
the  direction  opposite  to  the  positive  x-direction,  y  is  the 
mass  per  unit  length  of  the  rod,  and  L  is  the  length  of  the  rod; 

(*)  denotes  partial  derivative  with  respect  to  time,  and  (  )x  denotes 
partial  derivative  with  respect  to  position  x. 


To  answer  the  question  of  whether  or  not  there  exists  general 

variational  functions  corresponding  to  Equations  (1)  and  (2) ,  we 

[41 

can  use  the  simple  method  of  examining  the  integral  given  below  . 

Jl!  (3) 

If  it  is  possible  to  take  S  in  Equation  (3)  out  of  the  two  inte¬ 
grals,  then  it  means  that  the  general  variational  function  has  been 
found.  Otherwise,  it  does  not  exist.  Partial  integration  of 
Equation  (3)  gives 

- 4  Jl;  K  -H  <*>■-(-?-  <»->•-  -f ' 

+l'i'  ‘,i*  °  141 
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In  Equation  (4),  from  the  conditions  on  the  time  limits  of 
Hamilton's  principle,  the  third  term  vanishes.  The  fourth  term  also 
vanishes  because  of  Equation  (2)  and  the  fact  that  6cu  is  an 
allowed  function.  Therefore,  Equation  (4)  can  be  rewritten  as 

$L0^ir2*adxdt  (5) 

where 

Obviously,  H  is  the  general  Hamiltonian  of  the  system  when  g  is  a 
director  force  rather  than  a  follower  force.  For  the  Pfluger  rod, 
as  g  is  a  follower  force,  0  •  and  the  right  hand  side  of 

Equation  (5)  is  the  time-space  integral  of  the  virtual  work  done 
by  the  non-conservative  component  of  force,  -g&jx,  acting  through  the 
virtual  displacement  6<u.  It  can  thus  be  seen  that  the  general 
variational  function  cannot  be  obtained  for  the  system  from  ordinary 
energy  considerations.  Only  something  in  the  form  of  a  variational 
statement,  as  exemplified  by  Equation  (5)  ,  can  be  derived.  However, 
it  can  readily  be  seen  from  the  fact  that  the  right  hand  side  of 
Equation  (5)  is  hyperbolic  in  cux  and  5co  and  not  an  ordinary  second- 
order  expression,  that  the  expansion  of  the  approximate  solution 
will  always  lead  to  the  appearance  of  asymmetrical  matrices. 

Derivation  of  Equation  (5)  shows  us  that  if  we  can  somehow  take 
the  term  in  Equation  (1)  containing  out  of  the  process  of  partial 
integration,  then  it  will  be  possible  to  obtain  a  general  variational 
function.  In  fact,  substituting  5o)  with  6co  in  Equation  (3),  we 
obtain 

which  becomes,  after  partial  integration, 

rf*  r  1  -  *6 + -£-  a  -  *  >  d*dt 


(6) 


(7) 


1  J!|  5L.  -H  <-J‘-  u-  «)<»»)■ 

~J*  *»*■  l‘ “"+ Jt  (■ -S-  » 


Following  the  same  reasoning  as  above,  we  know  that  the  last  three 
terms  of  Equation  (7)  should  vanish  because  of  the  conditions  on  the 
time  limits  and  the  boundary  condition  given  in  Equation  (2).  Fur¬ 
thermore,  for  the  sake  of  convenience,  we  multiply  both  sides  of 
Equation  (7)  by  2,  and  obtain  the  generalized  Hamiltonian  principle: 

bH*=  0  (8) 

The  expression  for  the  generalized  Hamiltonian  variational  function  is 

{(WJ* — o) 

I 

Thus  we  have  established  the  variational  principle  corresponding 
to  the  given  differential  Equation  (1)  and  its  boundary  condition 
(2) .  The  approximate  solution  can  be  expanded  on  the  basis  of 
Equations  (8)  and  (9);  particularly,  the  finite  element  approximate 
solution.  Compared  with  an  expansion  based  on  the  variational 
statement  (5),  that  based  on  the  generalized  Hamilton's  principle 
makes  the  numerical  computation  easier  because  the  symmetry  property 
is  retained  in  the  latter. 

By  obtaining  the  generalized  Hamiltonian  variational  function 
and  the  generalized  Hamilton's  principle  we  have  in  effect  trans¬ 
formed  the  problem  of  the  non-conservative  system  into  that  of  a 
generalized  conservative  system.  Here  the  "energy"  cannot  possess 
any  physical  meaning  but  is  only  a  mathematical  analog.  In 
fact.  Equation  (3)  is  apparently  an  expression  of  the  principle 
of  virtual  work.  On  the  other  hand,  in  Equation  (6) ,  the  force  is 
multiplied  by  a  term  that  is  not  its  corresponding  virtual  displacement ,- 
and  the  product  is  not  virtual  work  in  the  physical  sense.  This  type 
of  mathematical  virtual  work  or  energy  can  be  termed  generalized 
virtual  work  or  energy.  Compare  this  with  the  canonical  problem  of 
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conservative  systems.  In  the  latter,  the  general  Hamiltonian 
variational  function  is  the  time- integral  of  the  Lagrangian  operator 
which  is  the  difference  between  the  kinetic  energy  and  the  potential 
energy  of  the  system,  and  the  total  energy  of  the  system  (sum  of 
the  kinetic  and  potential  energy)  is  conserved  (i.e.,  its  rate  of 
change  with  time  is  zero) .  Similary,  from  the  generalized 
Hamiltonian  variational  function.  Equation  (9),  we  can  write  (twice) 
the  generalized  total  energy  of  the  system  as  follows: 

do) 

Differentiate  Equation  (10)  with  respect  with  time  to  obtain 

M- -  2  |  w  +  — (£.—  •*)  w^-w.Jx 

+  2(WW.+-J-U,„w„)|^  (11) 

From  Equations  (1)  and  (2),  M=0,  i.e.,  with  respect  to  the  generalized 
total  energy  M,  the  system  is  constant,  steady  and  conservative. 
Therefore,  we  may  say  that  the  original  problem  of  the  non-conservative 
system  has  been  transformed  into  that  of  a  generalized  conservative 
system. 


Let 

w(x,  *)  d2) 

Substitute  this  into  Equations  (8)  and  (9),  and  making  use  of  the 
conditions  on  the  time  limits,  we  obtain  the  following  boundary  value 
variational  expression; 

&A-  0 

(13) 

where  / 


y: 


a 


;s 

l\V 


a 


N  , 

s  • 


I 


►  *  • 


u 


Equations  (13)  and  (14)  are  the  starting  points  for  establishing 
the  finite  element  model.  Similarly,  the  boundary  condition  (2) 
is  changed  to 


»r(0)-»'a)-»r-(0)-W'.U)-  0 


The  differential  Equation  (1)  will  take  the  form 


— a>W H — -  x)W„=  0 


(15) 


(16) 


Based  on  the  generalized  Hamilton's  principle,  we  may  let  co=0. 
and  determine  the  critical  load  by  applying  the  conventional  iterative 
technique  for  Rayleigh  quotient.  In  order  to  illustrate  this, 
we  will  briefly  discuss  the  form  of  instability  of  Pfluger's  rod. 
Consider  the  following  general  function 


F( ®*.  9,rr(x))  =  +  —~(L  -  x)WmSjWw4x=  o  (1 


7) 


From 


dF d9+^riW- 0 

am1 


where 


(18) 


it  can  be  checked  by  partial  integration  that  the  two  terms  on  the 
right  hand  side  of  Equation  (18)  are  equal.  Furthermore,  from 
Equation  (16)  we  know  that  aF/aW  *dW—  0  .  Therefore,  the  slope 

of  the  eigenvalue  curve  . 


dca1  aF  lag _ J  o  ^ 


1 


L  U-x ) 


{Waydx 


dg 


-WW„dx 


-> 
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From  Equation  (15),  we  know  that  the  denominator  of  the  above 
expression  is 


Therefore, 


JLt  ovydx 


dto1 


(19) 


dg 


<  0 


Equation  (19)  shows  that  the  slope  of  the  eigenvalue  curve  of  Pfluger' 
rod  is  always  negative.  Therefore,  instability  will  occur  in  the 
form  of  divergence  rather  than  flutter.  Hence,  if  one  is  only  con¬ 
cerned  with  the  determination  of  the  critical  load,  one  may  let 
co=  0,  and  the  problem  becomes  one  of  stability  of  static  forces. 

It  will  not  be  hard  to  see  that  the  problem  of  determining  the 
critical  load  g  . .  for  the  function  W(x)  that  satisfies  the 
boundary  condition  (15)  is  equivalent  to  that  of  determining  the 
minimum  value  of  the  following  Rayleigh  quotient. 


, - N(JV')_ 

Dim 


r. 


(tv„ydx 


(20) 


Obviously,  the  numerator  and  denominator  of  the  above  expression 
are  positive  definite  integrals.  The  value  of  G  determined  from  the 
stationary-value  condition  on  Equation  (20)  is  nothing  but 
Equation  (13)  for  oj=  0. 


III.  FINITE  ELEMENT  MODEL  FOR  PFLUGER  ROD 


Compared  with  the  finite  element  discretization  based  on  a 
variational  statement,  that  based  on  the  generalized  Hamilton's  varia¬ 
tional  principle  formulated  above  has  the  following  characteristics 
besides  possessing  symmetry.  First,  the  finite  element  discretiza¬ 
tion  model  has  to  satisfy  all  boundary  conditions  and  not  just  the 
geometric  boundary  conditions  of  the  original  problem.  Secondly, 
the  requirement  on  the  order  of  continuity  among  the  elements  is 
increased.  The  highest  order  derivative  in  the  general  variational 
function  is  of  the  third  order  while  those  for  the  method  using 


variational  statements  and  for  the  usual  problem  of  conservative 
systems  are  only  of  the  second  order. 

With  respect  to  the  second  characteristic,  it  would  seem  that 
the  so-called  "over-conservative  connected  elements"  would  be  useful 


in  problems  of  stability  of  non-conservative  systems.  However, 
available  over-conservative  connected  elements  are  not  adequate  for 
forming  c  omplete  element  models  because  up  until  now  there  is  as 
yet  no  element  used  in  structural  analysis  that  has  been  derived  from 
a  general  function  containing  third  order  derivatives.  Moreoever, 
the  existence  of  distributed  follower  forces  results  in  complicated 
generalized  geometric  stiffness  matrix. 

For  the  sake  of  convenience,  we  express  Equation  (14)  in 
dimensionless  quantities.  Let 


Then,  Equations  (14)  and  (15)  can  be  rewritten  as 

S"  J  3  {QW;)*- (»';:;)*+  9  (  1  -  £  )  {IV, .y ')  <*C 
(  0  )  *  (  1  )  =  (  0  )  -  (  1  )  -  0 


Let  the  rod  be  divided  into  n  elements  of  equal  length.  The 
elements  and  nodes  are  numbered  from  1  to  n  and  from  1  to  (n+1) , 
respectively,  along  the  positive  x-direction.  Let  8  denote  the 
local  coordinate  of  each  element,  with  values  lying  between  0  and  1. 
After  carrying  out  the  finite  element  discretization,  the  general 
function  h  will  be  in  the  form  of  an  algebraic  expression.  After 
completing  the  variational  operations,  we  obtain  the  algebraic 
eigenvalue  problem  to  be  solved: 

tK<-9{Kt-Ks)UD}=QKKj{D}  (24) 
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where 


{£)}  —  is  the  effective  degree  of  freedom 

—  is  the  generalized  mass  matrix 

I 

»  V  ( (»  —  »  +  1 ) —  is  t*le  generalized  geometric  stiffness  matrix 
i  (corresponding  to  the  conservative  force) 


C KJ=  n  ^CATS) 


—  is  the  generalized  geometric  stiffness  matrix 
(corresponding  to  the  non-conservative  force) 


C Kfi  —  is  the  generalized  stiffness  matrix 

I 

iKfi,  1,  2,  3,  4,  —  denotes  the  matrix  related  with  the  ith  element. 


This  matrix  will  have  different  forms  of  expression  depending  on  the 
different  design  elements  of  the  shape  function.  However,  the  form 
will  not  vary  with  the  element  number  i.  Therefore,  one  may  omit  the 
superscript  i.  The  set  of  matrices  listed  below  are  the  shape 
functions  corresponding  to  fifth-order  complete  polynomials.  The 
nodal  variable  of  the  ith  element  is 


uyj==L  OV,  W„)„  Or,  Wu  w„)Ul  j 


Even  though  if  one  is  only  concerned  with  the  critical  load,  is  of 
no  use,  we  list  its  values  below  for  the  sake  of  completeness  and 
for  later  application  in  the  problem  of  the  rectangular  plate. 
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270 

288 

symmetrical 

rH  ' 

II 

* 

15 

21 

2 

1260 

-1800 

-270 

-15 

1800 

270 

-18 

-6 

-270 

288 

>  “15 

6 

1 

15 

-21 

2  j 

'  3600 

'I 

1800 

1152 

i 

1 

KI--1 

90 

66 

18 

210 

-3600 

-1800 

-90 

3600 

1800 

648 

24 

-1800 

1152 

>  -90 

-24 

3 

90 

-66 

18  j 

'  1200 

'l 

480 

264 

jh-4 

40 

24 

3 

140 

-1200 

-480 

-40 

1200 

720 

216 

16 

-720 

504 

-20 

0 

1 

20 

.  -20 

9  / 

f 

2880 

1440 

768 

240 

144 

36 

4  4 

- 

2880  - 

1440 

-240 

2880 

1440 

672 

96 

-1440 

768 

-240 

-96 

-12 

240 

-144 
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IV.  FINITE  STRIP  MODEL  FOR  SIMPLY  SUPPORTED  RECTANGULAR  PLATE 


Let  the  boundaries  of  the  simply  supported  rectangular  plate  be 
at  x=0,  a;  and  y=0,b.  The  plate  is  under  the  action  of  tangentially 
distributed  follower  forces  along  the  negative  x-direction.  Its 
differential  equation  and  boundary  conditions  are: 

D.E.:  w  +  ^-VV+-f-  («-*)»_-  0  (25) 


B  .C.  : 


*  “  0  ,  a  ,  w  —  o  ^in+wn«  o 

^  “  0  ,  6  |  “  0  \TCm’m  o 


(26) 


Obviously,  there  exists  no  general  Hamiltonian  variational 
function  in  the  ordinary  sense  for  this  problem  because  the  external 
force  does  not  do  actual  work.  Following  the  method  used  for 
Pfluger's  rod,  we  multiply  the  right  hand  side  of  the  differential 
Equation  (25)  by  <$<o  and  perform  integration  over  space  and  time. 
Noting  that  not  only  is  the  bending  moment  zero  but  also  m  =0, 

XX 

and  on  the  perimeter  of  the  simply  supported  plate,  we  can 

derive  the  following  generalized  Hamiltonian  variational  function  and 

variational  principle:  . 

OH  *■  ® 

J'.. (27) 

■*  f'C®  —  Jf  )(»«)*  jdxdydt  ^8) 

Similarly,  we  can  also  demonstrate  the  conservative  property  of  its 
generalized  total  energy,  its  instability  being  in  the  form  of 
divergence,  the  existence  of  a  generalized  Rayleigh  quotient,  etc. 
These  will  be  omitted  here. 


y ,  t)-*'~rr(x,  y) 

from  Equation  (27),  Equation  (28)  can  be  rewritten  as 

&A-  0 

J!  uwj)'+ujv„y+(w^ 

+-±-(a-x)(jrj']dxJy 


The  boundary  condition  changes  to  the  following: 


On  the  perimeter. 


Let  *  y  ~65i  ,  X  be  the  ratio  of  the  length  of  the  plate 

to  its  width  such  that  b  =  Xa,  and 


Q1  =  (l<D*a«/.D 
=9a*/D 


(33) 


Then  Equations  (31)  and  (32)  are  changed  correspondingly  into  their 
dimensionless  forms: 

i-  f1  Jl  {Q*(lK0*+[(^;K)«+-^-(^«J*+-^r  (^ui)*] 

(34) 

On  the  perimeter,  where  £  =0.1  and  5  =  0.1,  we  have 

0  (35) 

The  finite  element  model  for  the  rectangular  plate  can  be 
expanded  on  the  basis  of  Equations  (34)  and  (35).  Here,  we  adopt 
the  convenient  finite  strip  model,  and  make  use  of  the  matrices 
that  are  already  available  for  the  element  model  of  Pfluger's  rod. 


Let 

x,y)=m  ^  F.(  £  )sinmjt| 

■t »  » .  j  .  s  — 

(36) 


Substitute  Equation  (36)  into  Equation  (34),  and  complete  integration 
with  respect  to  Note  the  orthogonality  of  the  trigonometric 

functions.  It  follows  that  for  different  values  of  m,  there  is  no 
coupling  among  the  corresponding  F  (s)'s.  We  obtain 


+ ?  ( i  —  'j; 


(37) 


For  the  discretization  along  the  ^direction,  we  adopt  a  fifth  order 
polynomial  as  in  the  case  of  Pfluger's  rod,  and  divide  the  plate 
into  n  equal  elements.  Comparing  Equations  (37)  and  (22),  we  see 
that  for  any  given  odd  value  of  m,  the  Fm  in  Equation  (37)  corres¬ 
ponds  to  the  W  in  Equation  (22),  except  that  the  derivative  of  Fm 
with  respect  to  ?  is  denoted  by  and  the  expression  corresponding 

to  (W^  )  2  is  made 

Hence,  corresponding 
structural  matrices: 


kd  WW  V-*  ***2  f  k4*iV4  V  ny  i,  S  ij  i  VU  W4. 

up  of  three  terms  [  (F:)H  (FI)1  (ms)*  +  -L. ] 

ig  to  each  value  of  m,  there  are  the  following  J 


V”  c  o  V’ 


i 

CA'.JsCA'0  -  ■  2  ( «  +  1  -  O  CATJD 


-  is  the  generalized  mass  matrix 

-  is  the  generalized  geometric 
stiffness  matrix  (corresponding  to 
the  conservative  force) 


CATj.3  — CATO*  »  2  CATJ)  _  the  generalized  geometric  stiffness 

I 

matrix  (corresponding  to  the  non-conserva¬ 
tive  force) 

CKtm)  =  C  Kt  3  +  2  C  KJ  +  (—)*  C  KJ 

(38) 

=«‘2  cATi:+2  (-^Ljn2(n  +  i-ocAi:+(^V 


In  the  above  AT{,  ATi,  AT$,  Ki  are  the  same  as  those  for  Pfluger’s  rod 
with  the  same  notations. 

Finally,  we  present  some  simple  numerical  examples.  For  Pfluger's 

rod.  Reference  [5],  gives  qcr^t=18 .96 .  From  our  finitel  element  model, 

q  ■ ..=19 .188  for  n=l  and  q  =1^97  for  n=2 .  For  the  simply  supported 
cn  l  cn  u 

plate,  taking  A=n=m=l,  we  find,  from  Equation  (38)  qcrit=68 . 16 .  The 
approximate  solution  given  in  Reference  [2]  is  69.1. 


V.  CONCLUSION 

We  have  used  simple  approach  to  discuss  Pfluaer's  rod  in  detail, 
particularly  concerning  the  variational  principle  that  is  used  to  trans¬ 
form  this  non-conservative  problem  into  one  of  a  generalized  con¬ 
servative  system.  On  the  basis  of  this  variational  principle,  we 
have  carried  out  finite  element  discretization  tests  and  presented 
the  element  matrices  that  are  all  symmetrical  and  convenient  to  use. 

We  have  also  pointed  out  the  procedure  by  which  the  finite  strip 
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model  for  the  related  simply  supported  rectangular  plate  may  be 
obtained  using  these  matrices.  It  should  be  pointed  out,  however, 
that  it  may  not  be  possible  to  transform  a  more  general  non-conserva¬ 
tive  problem  in  the  same  manner,  and  that  this  more  general  system 
may  have  some  new  characteristics  that  need  to  be  dealt  with. 

All  these  problems  require  further  study. 
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BUZZ  IN  AXISYMMETRIC  SUPERSONIC  INLET  AND  ITS  CONTROL 


Ho  Chung-wei  and  Chang  Shih-ying 
(Nanching  Aeronautical  Institute) 


ABSTRACT 


An  investigation  on  buzz  in  an  axisymmetric  external  compression  inlet 
and  its  control  are  presented.  A  zone  of  boundary  layer  bleed  on  the  cone 
surface  at  throat  has  been  designed  (Fig.l).A  block  diagram  of  a  dynamic 
pressure  measuring  system  employed  is  shown  in  Fig. 2. The  static  pressure 
fluctuation  on  the  cone  surface  at  throat  has  been  picked  up  by  Kulite  tra¬ 
nsducer.  The  results  of  the  experimental  investigation  can  be  summarized  as 
followst 

1.  At  Mach  number  and  angle  of  attack  <*“0",  there  is  only  a 

little  static  pressure  fluctuation  in  the  stable  supercritical  regime  or  subcritical 
regime  with  a  bow  shock  attached  ahead  the  cowl  lip.  When  the  bow  shock  is 
away  from  the  cowl  lip,  the  feature  of  anstability  appears.  And  as  buzz  occurs 
in  the  inlet,  the  pressure  fluctuates  sharply  (Fig. 3). The  characteristic  frequ¬ 
ency  of  pressure  fluctuation  is  43.75Hz  (Fig. 4). Its  A pt.mu.  is  about  several 
times  that  of  the  supercritical  regime  (Fig.  4). 

2.  At  M-“l. 97, and  <*“12*,  as  buzz  occurs  in  the  inlet,  the  bow  shock  on 
the  lee  side  of  the  cone  vibrates  sharply  with  a  large  amplitude  (Fig. 6),  and 
the  pressure  fluctuation  is  move.  sharp’'  than  that  on  the  windward  side 
(Fig.  7),  and  its  main  frequency  is  18.75Hz  (Fig. 8). The  Ap,...*  is  about  dou¬ 
ble  that  on  the  windward  side.  The  bow  shock  on  the  windward  side  vibrates 
with  a  small  amplitude  (Fig.S),  and  the  main  frequency  is  20Hz. 

3.  By  opening  bleeder  door  with  mass  flow  ratio  W h/W ,“0.0061,  the  inlet 
buzz  disappears  at  once.  The  bow  shock  becomes  stable  (Fig.  10),  and  the 
amplitude  of  the  pressure  fluctuation  and  their  A become  smaller  (Fig- 
10  and  Fig.  11). The  main  frequency  of  pressure  fluctuation,  for  example, 
on  the  lee  side  disappears  (Fig.  12), 


NOTATIONS 


M 

Moo  Mach  number  of  on-coming  stream 

Ms  Mach  number  of  stream  on  surface  of  cone 

?Too  Total  pressure  of  on-coming  stream 

W  Inlet  mass  flow  of  capture 

Ac  Inlet  area  of  capture 

Aj_  Area  of  inlet 

<j>  =Aoo/Aj_,  flow  parameter 

x  Distance  on  the  cone  from  the  point  of  measurement  to  the  apex 
PA  Density  of  power  spectrum 
a  Angle  of  attack 

Ps  Static  pressure  of  stream  on  surface  of  cone 

Ap  Root-mean-square  value  of  fluctuating  pressure 

i  f  m  f  s 

Amount  removed  by  boundary  layer 

Aoo  Free  stream  area 

Rc  Radius  of  capture 

x  Distance  from  point  of  measurement  to  cowl  lip 
oj  Frequency 

Subscript  "sup"  denotes  stable  supercritical  condition 

I .  FOREWORD 

Supersonic  air  intake  buzz  has  always  been  a  topic  of  great 
interest  in  other  countries  ^  Research  reports  have  been 

presented  on  the  characteristics  of  inlet  buzz,  development  of  buzz, 
and  inlet  and  throat  flow  chracteristics  under  inlet  buzz  condition 
when  the  angle  of  attack  is  0°.  However,  aircraft  flight  tests  have 
shown  that  inlet  buzz  usually  occurs  under  flight  conditions  where  the 
angle  of  attack  is  relatively  large.  Therefore,  the  study  of  this 
problem  is  of  practical  significance. 


II.  EXPERIMENTAL  MODEL 


The  model  is  a  double-cone  axisymmetric  supersonic  inlet  con¬ 
structed  on  the  scale  of  1:9.94  (Figure  1).  The  angles  of  the  cones 
are,  respectively,  35°  and  50°.  The  angles  of  the  inner  and  outer 
lips  of  the  cowl  are,  respectively,  10°  and  14°.  The  radius  of 
cature  Rc  of  the  inlet  is  47.4  mm.  The  geometric  throat  of  the 
inlet  is  located  at  a  position  where  x/Rc=0.232.  The  conic  body 
of  the  model  is  hollow.  The  boundary  layer  that  has  been  removed 
from  the  surface  of  the  cone  is  sucked  into  the  hollow  inside  of  the 
cone  to  flow  past  the  branch  rod  and  flow  meter,  and  then  returned  to 
the  low-pressure  cabin  of  the  wind  tunnel.  A  series  of  static 
pressure  holes  have  been  designed  on  the  upper  and  lower  surfaces 
of  the  conic  body.  At  the  position  where  x/R  =8.842  is  placed  a 
rotating  rake  with  12  points  for  measuring  total  pressure.  The 
position  of  the  rotating  rake  during  tests  is  monitored  by  means  of 
a  special  indicator.  The  inlet  mass  flow  is  controlled  through 
the  tail-plug  cone. 


Key:  1-Throat;  2-  Rotating  rake;  3-  Tail  cone; 
of  Kulite  transducer;  5-  Unit:  mm 


4-  Position 


III.  WIND  TUNNEL  AND  MEASURING  SYSTEM 

The  tests  were  carried  out  in  the  supersonic  inlet  wind  tunnel 
of  Nangching  Aeronautical  Institute.  The  dimensions  of  the  nozzle  are 
160  mm  x  160  mm,  and  its  maximum  angle  of  downward  tilt  is  15°. 

The  Reynold  number  Re  of  the  wind  tunnel,  that  has  as  its 
characteristic  dimension  the  inlet  radius  of  capture,  is 
2xl06  -  5  x  106  (for  Moo  =  1.6-2.63). 

The  steady-state  measuring  system  consists  of  SYD-1  transducer, 
XJ-100  detector  and  Model  LS-5  digital  recorder. 


j-  NJ-l 


•  D.C.4.3V 


y*  CF-500 
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Figure  2.  Block  diagram  of  dynamic  pressure  measuring 
system. 

1-  SF-72  data  amplifier;  2-  SR-50C  magnetic  tape 
recorder;  3-  NJ-l  recorder;  4-  SBE-20  two-trace 
oscilloscope;  5-  CF-500  power  spectrum  analyzer; 

6-  CX-446  x-6  recorder;  7-  FV-624  filter;  8-  R.M.S. 
voltmeter;  9-  PZ-8  digital  voltmeter. 


The  dynamic  pressure  measuring  system  is  shown  in  Figure  2, 

NO .CQL-030-100  and  NO .CQL-070-25  Kulite  transducers  are  used  to 
measure  the  fluctuating  pressure  at  the  geometric  throat  and  on  the 
surface  of  the  cone.  The  electric  signals  put  out  by  the  trans¬ 
ducers  pass  through  the  SF-72  amplifier,  and  are  recorded  by  the 
SR- 50  C  magnetic  tape  recorder.  At  the  same  time,  these  signals 
are  monitored  by  means  of  the  SBE-20  oscilloscope.  Or,  the  NJ-1 
recorder  is  used  to  record  the  time  history  of  the  pressure.  The  /41 
signals  also  pass  through  the  FV-624  filter  and  R.M.S.  voltmeter, 
and  the  root-mean-square  value  of  the  fluctuating  pressure  is 
obtained.  To  obtain  the  principal  characteristic  frequency,  use 
CF-500  to  analyze  the  power  spectrum. 

IV.  DESIGN  OF  BOUNDARY  LAYER  BLEED  ZONE 

According  to  References  [2]  and  [4],  upon  separation  of  the 
boundary  layer  on  the  (slanted)  cone  surface,  the  inlet  flow  is 
hindered,  and  the  unstable  flow  of  the  boundary  layer  gives  rise 
to  the  inlet  buzz.  For  M  00  =1.97  and  a=0°,  we  used  the  procedure  given 
in  Reference  [6]  to  compute  the  development  of  the  boundary  layer 
along  the  surface  of  the  cone.  Our  results  show  that  the  boundary 
layer  develops  very  fast  on  the  surface  of  the  cone  near  the  throat, 
increasing  from  0.677  to  over  1.088  mm.  Furthermore,  the  Mach 
number  of  the  surface  of  the  cone  near  the  throat  obtained  under 
similar  conditions  under  supercritical  conditions  reaches  above 
1.40-1.77,  and  the  boundary  layer  of  the  second  compression  surface 
becomes  disordered  (5] , [6] .  Hence  when  bow  shock  appears  in  front 
of  the  lip  of  the  inlet,  the  intensity  of  the  shock  wave  is  sufficient 
to  cause  the  boundary  layer  at  that  location  on  the  cone  surface  to 
separate  [7].  To  control  inlet  buzz,  we  design  the  boundary  layer 
bleed  zone  on  the  surface  of  the  cone  near  the  throat.  Five  rows 
of  holes,  bored  vertically  into  the  surface,  are  arranged  alter¬ 
nately  along  the  perimeter.  The  diameter  of  the  holes  is  1.5  times 

the  thickness  of  the  local  boundary  layer.  The  total  area  of  the  holes 
2 

is  468  mm  ,  approximately  equal  to  1/10  of  the  area  of  the  geometric 
throat. 


Surface  for  0°  Angle  of  Attack 

For  14^=1.97,  ct=0  0 ,  the  NO . CQL-030-100  Kulite  transducer  measures 
fluctuation  of  the  static  pressure  on  the  cone  surface  at  the  geometric 
throat.  The  static  pressure  history  for  various  typical  mass  flow 
ratios  as  well  as  the  corresponding  series  of  shock  waves  under  steady 
supercritical  and  subcritical  conditions  are  shown  in  Figure  3.  It 
can  be  seen  from  Figure  3  that  under  steady  supercritical  and 
subcritical  conditions,  the  pressure  has  a  flat  time  history.  When 
the  bow  shock  is  farther  away  from  the  lip  of  the  inlet  (<p  =0.78), 


Figure  3.  Static  pressure  history  at  specific  mass 
flow  ratios. 


the  static  pressure  undergoes  marked  fluctuation  with  time.  When 
buzz  occurs  in  the  inlet,  the  fluctuation  becomes  very  pronounced, 
with  a  large  amplitude.  There  is  a  close  correspondence  between 


the  above-mentioned  characteristics  of  static  pressure  fluctuation 
and  the  measured  static  pressure  fluctuation  for  the 

various  mass  flow  ratios,  as  shown  in  Figure  4.  ($=0.78) .  The  value  o 
A p,.mj  is  several  times  that  of  the  steady  supercritical  regime 
when  buzz  occurs  in  the  inlet. 


Figure  4.  Apr.m.s  of  static  pressure  at  specific 
mass  flow  ratio 
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Figure  5.  Power  spectrum  of  static  pressure  on 

the  cone  surface  in  intake  buzz  conditions. 


\  v  “7 


Power  spectrum  analysis  using  CF-500  Power  Spectrum  Analyzer 
yields  for  the  time  history  of  the  pressure  under  intake  buzz  condi¬ 
tion  the  power  spectrum  PA-oj  for  the  fluctuating  pressure  under  / 

intake  buzz  condition. 


It  can  be  seen  that  the  principal  frequency  of  the  static  pressure 
fluctuation  under  intake  buzz  condition  is  43.75  Hz  (Figure  5). 


The  Dynamic  Characteristics  of  Inlet  Buzz  for  Larger  Angles 
of  Attack 


(1)  Characteristics  of  bow  shock  fluctuation  in  front  of  the 
inlet  lip 

Tests  on  the  dynamics  of  buzz  have  been  performed  under  the  con¬ 
ditions  moo  =1.97  and  a=12°.  By  slowly  moving  the  tail  cone,  one  can  se 
that  the  inlet  bow  shock  passes  buzz  characteristics  different  from 
those  for  0°  angle  of  attack.  Under  buzz  condition,  the  bow  shock  on 
the  lee  side  of  the  cone  undergoes  violent  vibrations  which  are 
invisible  to  the  naked  eye,  while  the  bow  shock  waves  on  the  windward 
side  are  clearly  visible  (Figure  6) .  Figure  6  has  been  obtained  by 
photographing  at  the  speed  of  1/50  seconds.  Owing  to  the  violent 
vibrations  of  the  bow  shock  on  the  lee  side  of  the  cone,  its  image 
on  the  Schlieren  photograph  is  blurred.  The  bow  shock  on  the 
windward  side, however,  yields  a  very  clear  image. 

(2)  Characteristics  of  static  pressure  fluctuation  on  the 
surface  of  the  cone  under  intake  buzz  condition. 

Figure  7  shows  the  leeward  and  windward  static  pressure  time 

history,  measured  by  means  of  N0.CAL-075-25  transducer  for  the  cone 

(at  the  throat),  that  corresponds  to  Figure  6.  It  is  clear  that, 

under  intake  buzz  condition,  the  leeward  static  pressure  fluctuation 

is  much  more  severe  than  that  on  the  windward  side.  Experimentally 

measured  fluctuating  pressure  Ap  und.er  the  same  conditions  is 

r,m,s 

twice  as  large  on  the  lee  side  as  on  the  windward  side. 
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(3)  Characteristic  frequencies  of  static  pressure  fluctuation 
on  the  surface  of  the  cone  under  intake  buzz  condition 


Spectral  analysis  of  the  pressure  history  on  the  leeward 
and  windward  sides  of  the  cone  by  means  of  the  CF-500  power  spectrum 
analyzer  yields  the  power  spectra  Figures  8  and  S  under  intake 
buzz  condition.  It  can  be  seen  from  these  that  the  principal 
characteristic  frequency  is  18.75  Hz  on  the  lee  side  and  20  Hz  on 
the  windward  side.  The  two  values  are  very  close  to  each  other. 

3 .  Improvement  of  the  Dynamic  Characteristics  of  the  Inlet  by 

Elimination  of  Intake  Buzz  Through  Boundary  Layer  Bleed 

(1)  Characteristics  of  inlet  bow  shock 

Buzz  suppression  tests  have  been  carried  out  under  the  condition 
that  Mco  =  1*97  and  a  =12°.  When  buzz  occurs  in  the  inlet,  the 
boundary  layer  bleeding  holes  are  opened,  and  both  the  leeward 
and  windward  bow  shocks  become  steady  (Figure  10) .  The  picture  was 
taken  at  a  speed  of  1/50  second.  The  large-amplitude  fluctuation  of  the 
bow  shock  on  the  lee  side  has  disappeared. 

(2)  Comparison  of  the  static  pressure  fluctuation  characteristics 
and  the  characteristic  frequencies  under  intake 

condition  with  those  after  elimination  of  buzz. 

After  buzz  is  eliminated,  there  are  significant  changes  in  the 

static  pressure  time  history  on  the  lee  and  windward  sides  of  the 

cone,  and  the  amplitude  of  the  fluctuation  decreses  dramatically 

(Figure  11) .  Experimental  measurements  show  that  the  static  pressure 

Ap  „  on  the  lee  side  after  buzz  is  eliminated  is  only  1/3  that  in 
r,m, s  J 

the  presence  of  buzz.  The  static  pressure  Apr  m  on  the  windward  side 
after  buzz  is  eliminated  is  only  1/2  of  that  in  the  presence  of  buzz. 
With  the  disappearance  of  buzz,  the  sharp  fluctuation  of  the  bow  shock 
on  the  lee  side  abates,  and  a  steady  condition  is  attained.  At  this 
moment,  the  characteristic  frequencies  of  the  static  pressure  fluctua¬ 
tion  shown  in  Figure  8  all  disappear  (Figure  12) . 


(M-  »  1.97  a  «12*  f  =0.9144  Wu/Wc  =  0 . C061) 


Figure  10. 


Schlieren  photograph  of  shock  waves  without  buzz. 


Figure  11.  Comparison  between  the  static  pressure  histories 
of  the  lee  side  and  windward  side  with  or 
without  buzz.  : 

Key:  1-  Lee  sidj;  2-  Windward  size;  3-  Buzz;  4-  Steady  state 
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4. 


Discussion 


The  fact  that  the  characteristics  of  inlet  buzz  for  larger  angles 
of  attack  are  different  from  those  for  0°  angle  of  attack  is 
closely  related  to  the  asymmetry  of  the  shock  waves  on  the  surface 
of  the  cone  at  larger  angles  of  attack.  This  asymmetry  leads  to  a 
mubh  higher  static  pressure  on  the  windward  side  than  on  the  lee 
side,  as  shown  in  Figure  13.  The  larger  the  angle  of  attack,  the 
greater  the  difference  in  static  pressure.  This  causes  the  boundary 
layer  on  the  windward  side  to  flow  toward  the  lee  side,  resulting 
in  an  increase  in  thickness  of  the  boundary  layer  on  the  lee  side. 

On  the  other  hand,  the  asymmetry  in  the  shock  waves  on  the  cone  sur¬ 
face  also  leads  to  a  much  higher  Mach  number  on  the  second  cone 
surface  on  the  lee  side  than  on  the  windward  side.  It  can  be  seen 
from  Figure  14  that  under  the  condition  that  Moo=19  7  and  a=  8°  or  12°, 
the  Mach  number  of  the  stream  on  the  cone  surface  on  the  lee  side  in 
front  of  the  cowl  lip  is  greater  than  1.6.  When  the  bow  shock  waves 
appear  in  front  of  the  cowl  lip#  their  intensity  is  already  high  /4 

enough  to  cause  the  thickened  boundary  layer  to  separate,  or  even 
enough  to  cause  separation  of  the  main  stream.  On  the  windward  side, 
however,  the  Mach  number  of  the  stream  on  the  surface  of  the  cone 
only  goes  up  to  about  1.25,  and  the  bow  shock  may  not  be  strong  enough 
to  cause  separation  of  the  boundary  layer.  This  difference  becomes 
more  pronounced  with  increasing  angle  of  attack.  This  explains  why  when 
buzz  occurs  in  a  symmetrical  inlet  at  a  large  angle  of  attack  the  bow 
shock  on  the  lee  side  goes  into  violent  large-amplitude  fluctuation. 

With  a  high  pressure  fluctuation  Ap  while  the  bow  shock  on  the 

r ,  m ,  s , 

lee  side  only  undergoes  small-amplitude  fluctuation,  with  low  pressure 

fluctuation  Ap„  „  . 

r  ,m,  s 

VI  CONCLUSION 

(1)  Our  study  on  the  dynamic  characteristics  of  the  static 
pressure  on  the  cone  surface  at  the  geometric  throat  of  a  double-cone 
axially  symmetric  supersonic  inlet  for  Moo  -1.9  7  and  a=0°  reveals 
that  when  the  inlet  goes  from  the  steady  supercritical  regime  into 
the  subcritical  regime,  and  finally  reaches  buzz  condition. 


the  static  pressure  time  history  gradually  changes  from  a  steady 

level  trace  into  one  with  obvious  static  pressure  fluctuation  which 

becomes  very  violent,  and  with  large  amplitudes,  under  buzz  condition. 

At  this  moment,  the  static  pressure  fluctuation  Ap  may  become 

xr  f  m  f  s 

several  times  that  under  supercritical  or  subcritical  conditions. 

The  principal  characteristic  frequency  of  the  static  pressure 
fluctuation  is  43  Hz. 

(2)  For  Mo,,  =1.9  7  and  a  =12°,  our  results  show  that  when  buzz 

occurs  in  the  inlet,  the  bow  shock  on  the  lee  side  of  the  cone  goes 

into  large-amplitude  fluctuation  first,  with  a  principal  characteristic 

frequency  of  pressure  fluctuation  of  18.75  Hz;  Ap  of  the  static 

it  f  m  f  s 

pressure  fluctuation  is  about  twice  that  on  the  windward  side.  The 
bow  shock  on  the  windward  side  only  undergoes  small-amplitude  fluctua¬ 
tion,  with  a  principal  characteristic  frequency  of  20  Hz,  close  to 
that  on  the  lee  side.  Once  the  boundary  layer  on  the  cone  surface 
(at  .the  throat)  is  removed  by  bleeding,  intake  buzz  disappears,  and 
there  is  marked  decrease  in  the  static  pressure  fluctuation,  with  the 

static  pressure  fluctuation  Ap  on  the  lee  side  equal  to  1/3,  and  th 

xr  f  m  f  s 

on  the  windward  side  equal  to  1/2  of  their  corresponding  values  under 
buzz  condition.  Furthermore,  as  the  buzz  is  eliminated,  the  violent 
large-amplitude  fluctuation  of  the  bow  shock  on  the  lee  side  disappears, 
and  the  principal  characteristic  frequencies  of  its  static  pressure 
fluctuation  all  disappear. 

(3)  Finally,  we  have  made  an  analysis  of  the  characteristics  of 
buzz  in  an  axi symmetric  inlet  at  a  relatively  large  angle  of  attack. 
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A  PREDICTION  METHOD  FOR  CHARACTERISTICS  OF  COOLED  TRANSONIC 
TURBINES  WITH  LOW  ASPECT  RATIO  AND  HIGH  LOAD* 

Chuang  Yii-nan,  Chao  Yu-ch'i,  Shih  Ching,  Han  Chien-yiien 

(Nanhua  Powerplant  Research  Institute) 

ABSTRACT 


At  present,  advanced  aircraft  tnrbines  make  a  feature  of  low  aspect 
ratio,  transonic  flow  in  cascade  passage,  high  aerodynamic  load  and  intensive 
cooling.  In  order  to  predict  their  characteristics  the  AMDC  formulaCl,u  for 
cascade  loss  is  corrected  on  the  basis  of  the  so-called  KQ  system  13]  as  in 
formulaC  3  ).  In  accordance  with  the  experimental  results  of  transonic  casca¬ 
des  obtained  in  our  country,  the  AMDC  formula  for  supercritical  wave  drag 
is  then  modified  into  formulas(  4  )and(  5  )for  transonic  and  high  subsonic 
turbines  respectively. 

For  the  sake  of  eliminating  the  cascade  choke  conditions  which  do  not 
exist  in  fact  but  often  occur  in  calculation,  a  suitable  relation  for  the  flow 
angle  at  cascade  exit  is  proposed. 

The  prediction  method  takes  into  account  the  various  operating  condi¬ 
tions,  including  subcritical, critical  and  supercritical  expansions,  limitfogsu- 
percritical  expansion  and  limiting  load,  and  a  computer  program  has  been 
written. 

The  characteristics  of  a  transonic  and  a  high  subsonic  turbine  model/ 
and  the  design  performances  of  a  cooled  turbine  with  low  aspect  ratio  have 
been  predicted  with  this  method,  and  compared  with  test  results  in  good 
agreement.lt  is  shown  that  the  accuracy  is  about  99  percent. 
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I.  INTRODUCTION 

Advanced  aircraft  turbines  are  designed  with  high  total 
temperature  of  burning  gas  inlet  and  high  load,  and  have  as  their 
characteristics  low  aspect  ratio,  high  load  transonic,  usually 
unshrouded  rotor  blades  and  intensive  cooling. 

We  have  introduced  the  KQ  system,  recently  announced  abroad, 
in  our  analysis  of  the  AMDC  formula  for  cascade  loss.  On  the  basis 
of  the  experimental  results  obtained  in  our  country  for  transonic 
plane  cascades,  we  have  modified  the  AMDC  formula  for  supersonic 
wave  drag,  and  have  compiled  a  complete  and  practical  set  of  empirical 
and  semi-empirical  formulas  for  loss  in  low-aspect  ratio,  high- 
load  transonic  cascades. 

To  solve  the  problem  of  the  appearance  of  false  choke  conditions 
when  the  AMD  [4]  cascade  exit  flow  angle  formula  is  applied,  we  have, 
the  basis  of  careful  studies  using  the  empirical,  semi-empirical  and 
approximate  theoretical  formulas  used  abroad  and  at  home,  chosen 

[5  ] 

the  Apohob  formula  that  is  more  suitable  for  high-load  transonic 
cascade  flow  characteristics. 

With  increasing  turbine  inlet  temperature,  the  relative  amount  of 
cold  flow  is  increased,  up  to  5-6%  for  small  turbines.  The  geometric 
throat  of  a  high-lead  transonic  turbine  is  usually  in  a  critical  or 
high  subsonic  state.  Therefore,  it  is  necessary  to  consider  the 
mixing  of  cold  flow  in  the  cascade  row  by  row.  A  formula  for  cooling 
flow  mixing  has  been  derived  from  the  theory  for  one-dimensional 
compressible  axial  flow. 

Finally,  a  DJS-6  source  program  for  computing  the  characteristics 
has  been  written  to  take  care  of  many  operating  conditions  that  may 
arise  such  as  subcritical  expansion,  mass  flow  choke,  supercritical 


expansion,  limiting  supercritical  expansion  and  limiting  load. 

II.  SYSTEM  OF  FORMULAS  FOR  CASCADE  LOSS  IN  TRANSONIC  TURBINES  WITH 
LOW  ASPECT  RATIO  AND  HIGH  LOAD 

Turbine  cascade  loss  is  usually  divided  into  blade  section  loss, 
secondary  flow  loss,  interstitial  loss,  wake  loss  and  super¬ 
critical  wave  drag  loss.  Within  the  past  30  years,  many  empirical  loss 
formulas  have  been  published.  The  most  famous  among  these  are  the 
Ainley/Mathieson/Dunham/Came  (AMDC)  system  of  loss  formulas.  / 

Table  1  gives  typical  results  obtained  from  calculations  using 
this  system  of  formulas 


TABLE  1.  Comparison  of  AMDC  predicted  efficiencies  with  test  results 
for  some  typical  turbines. 
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Key:  1-  Turbine  No.;  2-  Aspect  ratio  H/B;  3-  Stator  blade;  4-  Rotor 

blade;  5-  Mass  flow  coefficient  WA/u;  6-  Load  coefficient  AH/u2; 
7-  Pressure  ratio  tt^;  8-  AMDC  predicted  efficiency  n*; 

9-  Measured  ef f iciency  n %’>  10-Deviation  in  efficiency  An*; 

11-  Relative  deviation  or  test  efficiency  — . 


We  conclude  from  the  above  that: 

1.  For  ordinary  turbines,  there  is  good  agreement  between  the 
AMDC  predicted  efficiencies  and  the  measure  efficiencies.  This  is  the 
reason  why  this  system  of  formulas  was  widely  adopted  abroad  and  at  home 


in  the  70' s.  We  have  chosen  this  system  as  the  basis  for  our  system 
of  loss  formulas . 


2.  For  high  subsonic  turbines  with  low  aspect  ratio  and 
unshrouded  rotor  blades,  such  as  turbines  No.  2  and  No.  3,  the  AMDC 
predicted  efficiencies  are  lower  than  measured  values  by  about  2% . 

Here,  it  is  necessary  to  introduce  the  KQ  system. 

(1)  For  cascades  with  low  aspect  ratio,  the  AMDC  predicted 
secondary  flow  loss  is  on  the  high  side. 

For  ordinary  turbines,  secondary  flow  loss  increases  with 

decreasing  aspect  ratio.  This  relation  is  usually  generalized  as  an 
•  f  2 1 

inverse  relation  .  Based  on  this,  when  the  aspect  ratio  is  at 
a  low  value  of  0.383,  the  secondary  flow  loss  is  as  high  as  78.33%, 
higher  than  ordinary  statistical  values  ^ .  This  shows  that  the  relat 
should  be  modified  for  cascades  with  low  aspect  ratio  J. 


TABLE  2. 

Loss 

analysis 

for  turbine  No. 

2. 

l1'*  a  it 
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*  a  m  * 
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|  0.1236 
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0  0900  . 
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Key:  1-  Aspect  ratio;  2-  Secondary  flow  loss;  3-  Blade  section  loss; 

4-  Interstitial  loss;  5-  Total  loss;  6-  Stator  blade; 

7-  Rotor  blade. 

(2)  AMDC  predicted  interstitial  loss  for  unshrouded 
rotor  blades  needs  to  be  modified. 

Part  of  the  flow  does  not  participate  in  the  deflection  but 
leaks  through  the  radial  interstices  between  the  shrouded  rotor 
blades,  resulting  in  a  loss.  This  type  of  loss  is  also  present  in  un¬ 
shrouded  rotor  blades.  Furthermore,  for  the  unshrouded  rotor 
blade,  the  steep  pressure  gradient  between  blade  ba^in  and  the 
blade  back  produces  deflected  flow  at  the  blade  tip,  forming  a  tip 
interstice  vortex  which  meets  the  corner  vortex  there  rotating  in 


the  opposite  direction.  The  combined  flow  moves  downstream,  not  only 
affecting  the  low-speed  flow  in  the  boundary  layer  near  the  end  wall, 
but  also  disturbing  the  deflection  of  the  main  stream  in  neighboring 
regions.  (Refer  to  Figure  1.).  The  larger  the  relative  radial  inter¬ 
stices,  the  greater  their  effect  on  the  cascade  loss.  For  turbine 
No.  2  shown  in  Table  1,  when  the  radial  interstice  is  0,  the  rotor 
blade  total  pressure  has  a  recovery  coefficient  a  dropping  to  0.844. 

The  value  computed  from  the  AMDC  formula  for  interstitial  loss  is 
0.859.  This  needs  to  be  corrected. 

In  addition,  the  AMDC  formulas  do  not  include  the  effect  of  the 
compressibility  of  the  flow  on  the  loss,  nor  the  effect  of  the  extent 
of  choke  at  the  tail  edge  on  wake  loss.  Hence,  it  is  necessary  to 
introduce  the  KQ  system  to  correct  the  AMDC  .Loss  formulas. 

3.  The  efficiency  predicted  with  the  AMDC  formulas  for  transonic 
turbines  with  high  load  are  on  the  low  side.  The  formula  for  comput¬ 
ing  the  supercritical  wave  drag  needs  to  be  modified. 

To  increase  the  load  of  the  turbine,  the  Mach  number  at  the 
cascade  exit  of  high  transonic  turbines  is  generally  designed  at 
around  1.2.  Under  the  operating  condition  of  supercritical  expansion, 
a  series  of  interacting  waves  are  generated  inside  the  slanted  cut. 

When  strong  shock  waves  interfere  with  the  boundary  layer,  the  low- 
momentum  fluid  in  the  boundary  layer  passes  through  the  shock  waves  and 
slowed  down.  This  generates  a  reverse  pressure  gradiant  which  causes 
the  boundary  layer  to  separate,  and  the  loss  is  increased.  We  call 
this  the  supercritical  wave  drag  loss.  Under  this  condition,  the 
flow  spectrum  inside  the  slanted  cut  is  extremely  complex. 

In  the  AMDC  system,  (1) 

where  y^  is  blade  section  loss  coefficient,  y'  is  the  basic  blade 
section  loss  coefficient  when  the  friction  between  the  blade-surface 
and  the  end  wall  boundary  layer  is  taken  into  consideration,  sad  f(s) 
is  the  increase  in  loss  due  to  supercritical  wave  drag. 

/(*)-  1  +60(M,-  l)1  (2) 
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Figure  1.  Flow  field  on  free 
end  for  unshrouded  rotor  blade. 

Key:  1-  Relative  motion  of  the 
wall;  2-  Boundary  layer  flow; 
3-  Radial  interstice;  4-  Tip 
interstice  vortex;  5-  Corner  i 
vortex. 
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Figure  2.  Comparison  of  transonic 
cascade  test  results  with  AMDC 
prediction  [3] . 

Key:  1-  Formula  (4) ;  2-  AMDC  syst 

3-  Data  for  four  types  of  transonic 
cascade . 


Figure  2  gives  a  comparison  of  transonic  cascade  test  results  with 

[3] 

the  AMDC  predicted  values  .  Figure  7  gives  a  comparison  of  tests 
results  with  the  AMDC  predicted  values  for  a  high  subsonic  model 
turbine  operating  under  supercritical  expansion  condition.  This 
shows  that,  for  high  transonic  turbines,  the  AMDC  predicted  efficiencies 
are  on  the  low  side,  and  the  supercritical  wave  drag  formula  (2) 
needs  to  be  corrected. 

In  recent  years,  we  have  accumulated  in  our  country  quite  a  bit 
of  transonic  cascade  test  results  .  The  following  generali¬ 

zations  can  be  made  from  this  data : 

(1)  When  the  transonic  cascade  is  working  under  the 
designed  operating  condition,  although  the  exit  Mach  number  is  greater 
than  1,  in  general  around  1.2,  the  loss  is  not  very  high,  as  shown 
in  Figure  3.  Only  when  it  deviates  from  the  designed  operating 


.> 


condition  and  continues  to  undergo  supercritical  expansion 
is  the  loss  increased. 

(2)  The  variation  of  loss  with  Mach  number  before  the 
designed  operating  point  may  be  computed  from  the  KQ  system. 

This  conclusion  can  be  explained  by  comparing  Figures  4  and 
5.  At  the  transonic  design  point  not  only  do  the  shock 

waves  inside  the  slanted  cut  tend  towards  the  tail  edge,  but  also 
their  intensity  is  no  higher  than  that  of  the  shock  waves  for  M2=1.0. 
This  is  because,  when  the  transonic  cascade  operates  under  its 
designed  operating  condition,  the  supersonic  flow  inside  the  slanted 
cut  has  flow  characteristics  that  are  compatible  with  the  geometric  / 
properties  of  the  slanted  cut.  In  this  case,  the  blade  back  portion 
of  the  slanted  cut  is  not  a  source  of  disturbance  for  the  supersonic 
flow.  Although  the  shock  waves  at  the  tail  edge  become  stronger  with 
increasing  Mach  number,  at  the  design  point,  the  tail  shock  FG  inci-. 
dent  on  the  slanted  cut  and  its  reflected  wave  GI  are  weakened  by 
the  action  of  the  series  of  expansion  waves  CDE  and  their  reflected 
waves  in  the  slanted  cut,  (see  Figure  6) .  Under  the  designed  transonic 
operating  condition,  the  shock  waves  in  the  slanted  cut  are  thus  not 
as  strong  as  those  under  the  critical  condition.  In  this  case,  the 
loss  due  to  interference  of  the  boundary  layer  by  the  shock  waves 
is  not  as  that  under  the  critical  condition. 


Figure  3.  Variation  of  transonic  cascade  loss 
with  exit  Mach  number  [9]. 

Key:  1-  Design  point 


Figure  4 .  Schlieren  photograph 
of  M2  =  1.0  spectrum  for  a 
transonic  cascade. 

Key:  1-  KTB -stator  blade  root. 


Figure  6.  Flow  pattern  in  a 
transonic  cascade  at  M2=l . 2 . 


Figure  5.  Schlieren  photograph 
of  M2  *  1.2  spectrum  for  a 
transonic  cascade. 

Key:  1-  KTB-stator  blade  root. 
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Figure  7.  Comparison  of  pre¬ 
diction  with  test  results  for 
a  high  subsonic  model  turbine 
in  supercritical  expansion 
condition. 

Key:  1-  Experimental  curve; 

2-  Curve  obtained  from  compu¬ 
tation  using  Equation  (5) ; 

3-  Curve  obtained  from  computa¬ 
tion  using  the  AMDC  formula. 


When  the  basic  blade  section  loss  coefficient  is  computed 
from  the  formula  corrected  on  the  basis  of  the  KQ  system 


where  Yp^pc  is  the  basic  blade  section  loss  coefficient  of  the  AM  DC 
system,  Kp  is  the  correction  factor  to  account  for  the  compressibility 
of  the  flow,  and  y^ocfc  is  the  loss  coefficient  of  the  front-edge 
shock  waves,  then,  on  the  basis  of  the  AMDC  loss  formulas.  Equation  (2) 
is  modified  to  read 


f(s)  =  1  +  60  <M2-Maesign) 


In  the  above,  Mdeg^gn  is  the  design  Mach  number  of  the  transonic 
cascade. 

The  validity  of  this  formula  has  been  borne  out  by  tests  per¬ 
formed  on  our  high- load  transonic  model  turbine  (see  Figure  11) . 

For  subsonic  and  high  subsonic  turbines,  the  slanted  cut  under¬ 
goes  supercritical  expansion,  deviating  from  the  design  operating 
condition,  and  the  loss  due  to  supercritical  wave  drag  should  be 
taken  into  consideration.  However,  the  critical  Mach  number  at  the 
exit  of  an  actual  turbine  is  not  1.  We  modify  Equation  (2)  on  the 
basis  of  Equation  (3)  to  read 

/(•)-  1  +60(M,-Af„)*  (5) 

where  M  is  the  exit  critical  Mach  number, 
cr 

The  validity  of  the  above  equation  is  borne  out  by  tests  on 
total  performance  of  our  high  subsonic  model  turbine  (Figure  7) . 

See  Table  3  for  the  exit  critical  Mach  number. 


TABLE  3.  Exit  critical  Mach  number  for  a  high  subsonic  model 


turbine. 
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Key:  1-  Cascade;  2-  Parameter;  3-  Total  pressure  recovery  coefficient 

4-  Total  pressure  loss  coefficient;  5-  Tail-edge  deflection 
angle;  6-  Tail-edge  wedge  angle;  7-  Exit  structure  angle; 

8-  Exit  critical  Mach  number;  9-  Stator  blade; 10-  Rotor  blade 


III.  EXIT  FLOW  ANGLE  FORMULA 

All  the  formulas  for  the  exit  flow  angle  have  the  form  arcsin  v  t 
whether  they  are  approximate  theoretical,  empirical  or  semi-empirical 
formulas.  All  of  them  are  based  on  [4],  [5],  [11]. 
K  is  the  deflection  factor  determined  by 


(1)  exit  Mach  number, 

(2)  the  slanted-cut  geometrical  characteristics  determined 
by  the  geometric  exit  angle,  tail-edge  deflection  angle 
and  tail-edge  wedge  angle,  and 

(3)  loss  in  the  cascade  channel. 

We  performed  some  computations  on  the  turbine  characteristics^ 

[4] 

using  the  Ainley/Mathieson/Dunham  formulas 


W|*l  0i»™  aresta-^?- 

/ 

Mi<0.5  A»-l.I47«rc4n-7 - 11.1+ 4 

0.5<M»<1.0  4  (A,-0,JCl  -6M,+9M5-4MD 


The  false  choke  condition  shown  in  Table  4  and  Figure  8  was  generated 
in  the  calculations.  This  is  because  the  deflection  factor  used  in 
the  AMD  formulas  was  induced  from  subsonic  cascades  and  is  not  com¬ 
pletely  suitable  for  high  subsonic  near-critical  operating  conditions. 


-  .  * 
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Taking  into  consideration  the  high  transonic  cascade  flow  characteris¬ 
tics,  we  choose  to  use  the  approximate  theoretical  formula,  the  Apohob 
formula  that  was  derived  from  the  basic  equations  starting 

out  from  the  concept  of  passage  and  taking  loss  into  account: 


Pi./— arcsn^-^-^  P»— arcsin  (x-jir 


J 


K  —cos 


T  \  cos  p.  J  \  3,  J 

yo^pf) 


(6) 


TABLE  4.  Generation  of  false  choke  condition  with  AMD  cascade  exit 
flow  angle  formula 
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9.950ft 
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t 

9.  mu  j 

9.1799 

Key:  1-  Operating  condition;  2-  Turbine  No.;  3-  Critical; 
4-  False  choke. 


Figure  8.  Generation  of  false  choke  condition  with  AMD 
cascade  exit  flow  angle  formula. 

Key:  1-  Critical  mass  flow  curve;  2-  False  choke  condition; 

3-  Exit  flow  angle  curve;  4-  Critical  operating  condition 


IV.  COOLING  FLOW  AND  MAIN  FUEL  FLOW  MIXING  FORMULAS 


When  carrying  out  calculations  for  cooling  flow  mixing  in  the 
cascade  of  a  high-load  transonic  cooled  turbine  row  by  row  using  the 
one-dimensional  compressible  axisymmetric  flow  theory,  we  first 
calculate  the  velocity  triangle  and  aerodynamic  parameters  for  the 
exit  main  fuel  flow.  Then  cooling  flow  mixing  is  carried  out  in 
the  axial  interstices.  Finally,  the  entrance  parameters  for  the 
next  row  are  calculated. 

The  assumptions  made  in  the  derivation  of  the  cooling  flow 
mixing  formulas  are: 

(1)  Both  the  cooling  flow  and  the  main  burning  fuel 
flow  are  axisymmetric. 

(2)  The  cooling  flow  and  the  main  fuel  flow  are  mixed 
evenly  in  the  axial  interstices. 

(3)  Within  the  axial  interstices,  the  loss  due  to  viscous 
friction  within  the  flow  and  between  the  flow  and  the  end  wall  may  be 
neglected. 


(4)  The  flows  are  steady  constant  flows. 

The  cooling  flow  mixing  formulas  derived  from  the  one-dimensional 
compressible  axisymmetric  flow  theory  are  as  follows. 
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Table  5  gives  an.  example  of  computation  using  the  above 
formulas.  The  predicted  values  deviate  from  the  measured  values 
by  0.488%,  justifying  the  use  of  this  simplified  theory. 


TABLE  5.  Prediction  of  cooling  flow  mixing  for  the  turbine  No.  2. 
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Key:  1-  Parameters;  2-  Relative  mass  flow  of  cooling  flow; 

3-  Relative  change  in  total  temperature;  4-  Relative  change 
in  total  pressure;  5-  Relative  change  in  static  pressure; 

6-  Relative  change  in  velocity;  7-  Relative  change  in  axial 
velocity;  8-  Relative  change  in  Mach  number;  9-  Change  in 
flow  angle;  10-  Relative  change  in  static  temperature; 

11-  Calculated  efficiency  for  cooling  flow  mixing;  12-  Measured 
efficiency  for  cooling  flow  mixing;  13-  Calculated ; 


14- 


'measured' 


15-  Numerical  values. 


V.  COMPARISON  OF  PREDICTED  CHARACTERISTICS  WITH  MEASURED  RESULTS 


Table  6  and  Figures  9-11  give  a  comparison  between  the  calculated 
and  measured  characteristics  of  a  high  subsonic  model  turbine  and 
a  high-load  transonic  model  turbine  made  in  this  country.  The  error 
is  1%. 
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Table  5  gives  a  comparison  between  the  calculated  and  measured 
efficiency  for  a  cooled  high  subsonic  turbine  with  low  aspect  ratio 
at  its  design  operating  point.  The  error  is  0.448%.  / 

TABLE  6.  Characteristics  of  design  condition  for  transonic  and 
high  subsonic  model  turbines 
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I  _  t  «»•  :*•* 
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Key:  1-  Parameters;  2-  Expansion  ratio;  3-  Equivalent  work; 

4-  Stagnant  efficiency;  5-  Stator  cascade  exit  Mach  number; 

6-  Rotor  cascade  exit  relative  Mach  number;  7-  Flow  deflection 
angle;  8-  Rotor  blade  convergence;  9-  High  subsonic  model 
turbine;  10-  High-load  transonic  model  turbine. 
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Figure  9.  Comparison  of  predicted 
efficiencies  with  test  results  for 
a  high  subsonic  model  turbine. 


Figure  10.  Comparisons  of  pre¬ 
dicted  equivalent  flow  rate  with 
test  results  for  a  high  subsonic 
model  turbine. 

Key;  1-  Experimental  critical  poi 

2-  Calculated  critical  point; 

3-  Experimental  curve; 

4-  Calculated  curve 


Key:  1-  Experimental  curve; 

2-  Calculated  curve 


Figure  11.  Comparison  of  prediction  efficiency 
characteristics  with  test  results 
for  a  high-load  subsonic  model  turbine. 

Key:  1-  Critical  line  for  stator  blade;  2-  Critical 
line  for  rotor  blade;  3-  Notation  for  test 
points 

VI.  CONCLUSION  /5 

1.  We  have  shown  through  tests  on  the  total  characteristics  that 
the  modifications  we  made  on  the  AMDC  supercritical  wave  drag  formula 
are  justified.  The  modified  formula  is  applicable  to  high  subsonic 

or  transonic  turbines. 

2.  Better  fit  is  obtained  between  experimental  results  and  the 
characteristics  calculated  from  the  new  set  of  cascade  loss  formula 
and  exit  flow  angle  formulas .  This  shows  that  this  method  of  compu¬ 
tation  is  applicable  to  transonic  turbines  with  low  aspect  ratio 
and  high  load  as  well  as  to  ordinary  turbines. 

3.  It  has  been  shown  through  an  actual  computed  example  that 
it  is  appropriate  to  use  the  cooling  flow  mixing  formula  derived 
from  the  simplified  on^-dimens  onal  compressible  axisymmetric  flow 
theory.  However,  it  i  “in  to  be  studied  how  one  can  convert  the 
complicated  compound  cooling  flow  in  the  cascade  channel  into  the 
equivalent  exit  axisymmetric  flow. 


Notations 

Acr  Area  of  geometric  throat 
B  Chord  length 

$  Flow  angle 

$x  Angle  between  the  average  momentum  vector  at  the  throat 
and  the  bow  line 

cp  Specific  heat 

F  Turbine  passage  area 
G  Mass  flow 

H  Blade  height 

K  Deflection  factor 

k  Specific  heat  ratio 

e  Radius  of  curvature  of  blade  tail  edge 

M  Mach  number 

•  Throat  width 

AT*  Change  in  total  temperature  caused  by  cooling  flow  mixing 
Ap*  Change  in  total  pressure  caused  by  cooling  flow  mixing 
Ap  Change  in  static  pressure  caused  by  cooling  flow  mixing 
AW 2  Change  in  velocity  due  to  cooling  flow  mixing 

N  Ratio  of  rotational  speed  to  design  rotational  speed 
p  Static  pressure 

p*  Total  pressure 

R^  Mean  radius  of  axial  interstice 
T*  Total  temperature 

t  Static  temperature,  cascade  spacing 

W  Velocity  of  fuel  gas  flow 

x  Angle  between  normal  to  the  cross  section  of  the  throat  and 
the  mean  momentum  vector 

ip  Velocity  coefficient 

Yg  Ratio  of  the  velocity  component  of  the  cooling  flow  along  the 
perimeter  to  that  of  the  burning  fuel  flow 

Yx  Ratio  of  the  velocity  component  of  the  cooling  flow  along  the 
axial  direction  to  that  of  the  burning  fuel  flow 
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Notations  (Continued) 


AG  Mass  flow  of  the  cooling  flow 

AW^  Axial  velocity  change  due  to  cooling  flow  mixing 
AM2  Change  in  Mach  number  due  to  cooling  flow  mixing 
A&2  Change  in  exit  flow  angle  due  to  cooling  flow  mixing 
y  Total  pressure  loss  coefficient 
z,  Energy  loss  coefficient 
G  Ratio  of  mass  flow  to  design  mass  flow 


Subscripts : 

A  Axial 

1  Cascade  entrance 

2  Cascade  exit 

L  Cooling  air 

ef  Effective 

1-n  From  the  cascade  entrance  to  throat 
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Effect  of  Resin  Matrix  Toughness  on  Mechanical 
Properties  of  Composite  Materials 
Northwestern  Polytechnical  University 
Ma  Ruiran,  Zhou  Xizhen,  Wang  Yan,  and  Zhang  Wu 

ABSTRACT 

Using  the  diglycidyl  ether  of  biphenol  A— butadiene-acrylo¬ 
nitrile  rubber  with  random  carboxyl  groups — 2-ethyl-4-methylimi- 
dazol  epoxy  resin  system  as  a  representative,  the  effect  of  the 
resin  matrx  toughness  on  the  unidirectional  mechanical  proper¬ 
ties  of  the  glass  fiber  composite  material  was  investigated. 

The  experimental  results  showed  that  with  increasing  resin  matrix 
toughness ,  the  fracture  resistant  characteristics  of  the  compos¬ 
ite  material  obtained  increased  correspondingly.  However,  the 
amplitude  of  increase  is  smaller  than  the  increment  of  the 
toughness  of  the  resin  matrix  itself.  The  short  beam  shear 
strength  of  the  composite  material,  in  addition  to  the  tensile 
modulus,  and  the  transverse  mechanical  properties  also  improve 
with  increasing  resin  matrix  toughness.  The  effect  on  the  flex¬ 
ural  characteristics  and  longitudinal  mechanical  properties, 
however,  was  not  significant. 

INTRODUCTION 

In  recent  years,  the  use  of  rubber  to  toughen  epoxy  resin 

r  3/^3 1 

has  already  attracted  the  attention  of  many  people.  1  However 

the  rubber  used  was  mostly  a  butadiene  rubber  with  a  carboxyl 
group  at  the  end.  There  are  few  experimental  results  with  other 
types  of  rubber  reported.  Due  to  the  fact  that  butadiene  rubber 
with  a  random  carboxyl  group  has  already  been  produced  in  batches 
in  our  country,  we  chose  this  rubber  as  a  toughening  agent.  Its 
toughening  effect  on  epoxy  resin  and  the  effect  of  the  toughened 
epoxy  resin  matrix  material  on  the  properties  of  the  fiber  re¬ 
inforced  composite  material  were  studied.  The  toughening  effect 

of  butadiene  rubber  with  a  random  carboxyl  group  on  epoxy  resin 

f  4  51 

had  already  been  reported  by  us.  '  The  results  of  the  study 
showed  that  after  being  modified  by  butadiene  rubber  with  a 
random  carboxyl  group  and  biphenol  A,  the  fracture  surface  energy 
of  the  biphenol  A  epoxy  resin  was  increased  by  1~4  times.  The 


impact  strength  was  increased  by  2~10  times.  The  toughness  modu¬ 
lus  (the  area  under  the  tensile  stress-strain  curve)  was  increased 
by  0.5~11.5  times. 

In  this  paper,  the  effect  of  the  toughness  of  the  epoxy  resin 
matrix  on  the  mechanical  properties  of  a  unidirectional  glass 
fiber  composite  material  was  reported.  Furthermore,  the  method  of 
testing  the  fracture  toughness  of  a  composite  material  was 
briefly  introduced. 

II.  TESTING  OF  FRACTURE  TOUGHNESS  KQ 

Kq  is  the  opening  type  critical  stress  strength  factor  of  the 
composite  material.  We  used  a  double  edge  notched  tensile  frac¬ 
ture  method  in  our  determination.  The  dimensions  of  the  test 
specimen  are  shown  in  Figure  1. 


- 18'— 
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Figure  1 .  Double  edge  notched  tension,  specimen 

The  stress  strength  factor  K  is  calculated  according  to  the 
following  formula 

/C-Koi'T  (1) 

where  a-  the  stress  acting  on  the  two  ends  of  the  specimen 

a  -  the  effective  crack  length.  The  initial  crack  length 
adopted  in  this  experiment  was  aQ  »  3.5mm. 

received  in  September,  1982 

Y  -  the  correction  coefficient  for  the  geometric  shape  of  /61 
the  specimen,  which  is  corrected  according  to  the 


following  formula: 


K-l.M+0.36- 1±  — 2.I2(4l),+3.^2(-|l)•  (A) 

w  is  the  width  of  the  specimen. 

The  effective  crack  length  is  obtained  using  the  compliance 
calibration  curve  method.  Figure  2  is  the  dimensionless  compli¬ 
ance  calibration  curve  obtained  using  eight  pieces  of  specimens 

Oft 

with  different  — -C  values  at  two  thicknesses  of  1  and  2mm.  From 

w 

Figure  2  one  can  see  that  the  dimensionless  compliance  curves 
obtained  with  specimens  fabricated  using  1mm  and  2mm  thick 
plate  materials  could  coincide  each  other  very  well. 


Figure  2.  Compliance  calibration  curve 

1 .  specimen  thickness 

2.  x  1mm 

3.  o  2mm 

4.  C-compliance 

5.  E-tensile  modulus 

6.  B-thickness 


Figure  3.  Schematic  representation  of  determining 

1 .  instability  point 

2 .  crack  point 

According  to  the  curve  in  Figure  2,  the  P-5  curve  of  each 
formal  specimen,  and  equation  (1)  ,  it  is  possible  to  obtain  the 
value  of  K  corresponding  to  the  effective  crack  length  in  the 
crack  propagation  process.  From  these,  the  K-a  curve  is  prepared 
(Figure  3).  In  this  paper,  the  first  turning  point  on  the  curve 
is  the  crack  inception  point.  The  second  turning  point  is  the 
crack  instability  point.  The  value  K  corresponding  to  the  insta¬ 
bility  point  is  the  fracture  toughness  Kq. 

III.  RESULTS  AND  DISCUSSION 

A  unidirectional  fiber  plate  is  the  basic  composition  unit 
of  a  fiber  composite  material.  For  this  reason,  in  order  to 
understand  the  effect  of  the  toughness  of  the  resin  matrix  (using 
impact  strength,  toughness  modulus,  and  fracture  surface  energy 
as  the  evaluating  indicators)  on  the  mechanical  properties  of  a 
glass  fiber  composite  material,  three  types  of  resin-hardener 
systems  were  selected  as  the  objects  to  be  studied  in  this  work. 
The  corresponding  unidirectional  glass  fiber  composite  materials 
were  prepared.  The  fracture  resistant  properties  and  conventional 
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mechanical  characteristics  of  these  three  composite  materials  were 
analyzed  and  compared. 

The  compositions  and  formulations  of  these  three  resin-hard¬ 
ener  systems  are  shown  in  Table  1.  The  major  thermal-mechanical 
properties  of  cast  resins  are  listed  in  Table  2. 

The  matrix  material  A  was  the  unmodified  618  epoxy  resin/2- 
ethyl-4  methyl imidazol  system.  B  and  C  were  modified  by  a  random 
carboxyl  butadiene  rubber,  and  a  random  carboxyl  butadiene  rubber 
and  diphenol  A,  respectively.  By  comparing  the  data  in  Table  2, 
it  is  possible  to  see  that  the  fracture  surface  energy,  impact 
strength,  and  toughness  modulus  of  the  matrix  materials  A,  B,  and 
C  were  improving  gradually.  This  shows  that  the  ability  of  these 
three  materials  to  resist  crack  propagation  and  the  capability  of 
resisting  fracture  under  a  fast  and  slow  load  will  increase  ac¬ 
cording  to  the  sequence  mentioned  above.  The  thermal  distortion 
of  these  three  is  more  or  less  the  same.  Therefore,  it  can  be 
considered  that  the  matrix  materials  B  and  C  have  a  relatively 
higher  toughness.  The  toughness  of  C,  however,  is  higher  than 
that  of  B.  The  experimental  results  of  the  fracture  resistance  /62 
characteristics  and  conventional  mechanical  properties  of  uni¬ 
directional  glass  fiber  composite  materials  corresponding  to  the 
three  matrix  materials  are  described  in  the  following: 


ua*«ita  (HaSi  **«.  o.u~«.34a*/MQS) 
XIHTRiK  <X*a> 

«»A  (tfXRt) 

utexaiK) 


Table  1.  Formulation  of  resin-hardener  systems 

1.  raw  material 

2.  weight  ratio 

3 .  resin  matrix 

4.  618  epoxy  resin  (biphenol  A  type,  epoxy  value: 
0.48~0.54  equivalent/ 100  gram) 
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5 .  butadiene  rubber  with  a  random  carboxyl  group 
(industrial  grade) 

6.  diphenol  A  (chemical  reagent  grade) 

7.  2-ethyl-4-methylimidazol 


*f  *  a  a  a 
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Table  2.  Properties  of  cast  epoxy  resins 

1 .  testing  item 

2.  cast  epoxy  resin  e  2 

3.  fracture  surface  energy,  xlO  erg/cm 

4.  impact  strength,  kg.cm/cm2 

5.  toughness  modulus  kg.cm/cm3 

6.  tensile  strength,  kg/cm2  4  2 

7.  tensile  elastic  modulus,  xlO*kg/cnr 

8.  fracture  elongation  rate,  % 

9.  flexural  strength, kg /cm2  .  2 

10.  flexural  elastic  modulus  ,xl0  xg/cm 

-  11.  thermal  distortion  temperature,  °C 

12.  flexural  stress s 18 ,5kg/cm2 

1.  Fracture  Resistant  Characteristics  of  Unidirectional 
Glass  Fiber  Composite  Materials 

(1)  Fracture  Surface  Energy 

The  testing  results  of  the  fracture  surface  energy  of  a  crack 
propagating  along  the  fiber  axis  are  shown  in  Table  3. 
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A  i 

21.43 

4.21 

B  i 

15.32 

11.7 

C  ! 

27.33 

30.1 

Table  3.  Fracture  surface  of  unidirectional  E-glass  fiber 

composites 

1.  epoxy  resin  matrix 

2 .  composite  material 

3.  glue  content,  %  52 

4.  fracture  surface  energy,  xlO  erg/cm 

A  comparison  of  the  data  in  Table  3  shows  that  with  increas¬ 
ing  matrix  toughness  the  fracture  surface  energies  of  composite 
materials  B  and  C  increased  by  69%  and  270%  as  compared  to  that 
of  A,  respectively.  This  is  because  the  fibers  could  not  be 
aligned  straight  in  a  composite  material.  When  the  tip  of  a 
crack  passes  through  the  resin  matrix,  the  two  surfaces  near  the 
crack  zone  of  a  crack  are  still  connected  by  the  unbroken  fibers 
or  by  broken  fibers  which  have  not  been  pulled  out.  Hence,  the 
fracture  surface  energy  measured  by  a  cleavage  test  is  primarily 
consumed  in  the  splitting  of  the  matrix,  the  detachment  of  the 
fiber-matrix  interface,  and  the  fracture  of  the  fibers.  Therfore 
increasing  the  toughness  of  the  matrix  and  the  strength  of  the 
fiber-matrix  interface  will  help  the  fracture  surface  energy  of 
a  composite  material  to  increase. 

A  comparison  of  the  short  beam  shear  strength  of  the  three 
composite  materials  A,  B,  and  C  (see  Table  6)  shows  that  the 
bonding  force  of  the  fiber-matrix  interface  is  improved  after 


the  epoxy  resin  is  toughened  by  rubber  and  diphenol  A.  There¬ 
fore,  according  to  the  difference  in  the  matrix  toughness  and 
interface  strength  among  these  three  composite  materials,  it  is 
possible  to  predict  that  the  improvements  in  the  fracture  surface 
energy  are  the  results  of  improved  matrix  toughness  and  interface 
strength . 

A  further  comparison  of  the  fracture  energy  of  the  compos¬ 
ite  material  and  the  matrix  shows  that  with  increasing  matrix 
toughness,  the  amplitude  of  increase  in  the  fracture  surface 
energy  of  a  composite  material  is  smaller  than  the  increment  in 
the  matrix  itself.  This  may  be  attributed  to  the  thin  film  state 
of  the  resin  matrix  in  a  composite  material  since  the  toughness 
of  a  thin  film  resin  is  significantly  lower  than  that  of  a  block 
of  resin. ^ 

(2)  Impact  Strength 

The  impact  characteristics  of  the  three  composite  materials 
are  tabulated  in  Table  4. 
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Table  4.  Impact  strength  of  unidirectional  E-glass  fiber  com¬ 
posites 

1.  epoxy  resin  matrix 

2.  composite  material 

3.  glue  content,  %  2 

4.  transverse  impact  strength,  kg.cm/cm  ^ 

5.  longitudinal  impact  strength,  kg.cm/cm 

6.  mean  value 

7.  scatter  coefficient 


8.  mean  value 

9.  scatter  coefficient 

A  comparison  of  the  data  in  Table  4  shows  that  the  composite 
materials  have  a  relatively  higher  transverse  impact  strength 
(the  length  direction  of  the  transverse  specimen  is  perpendicular 
to  the  axial  direction  of  the  fiber)  after  they  are  toughened. 
Composite  materials  B  and  C  are  raised  by  39%  and  69%  as  compared 
with  that  of  A,  respectively.  The  increments  of  the  longitudinal 
strength  of  a  composite  material  (the  direction  of  the  length  of 
a  longitudinal  specimen  is  parallel  to  the  fiber  axis)  are  rela¬ 
tively  smaller.  B  and  C  improved  by  16%  and  23%,  respectively, 
as  compared  to  A. 

As  shown  in  Figure  4,  the  fracture  mode  of  a  composite  mate¬ 
rial  transverse  impact  specimen  is  primarily  the  detachment  of  the 
interface  and  the  cracking  of  the  matrix.  Therefore,  the  trans¬ 
verse  impact  strength  is  mainly  determined  by  the  tensile  proper¬ 
ties  and  interface  strength  of  the  matrix  material.  Because  the 
toughness  modulus  of  the  matrices  B  and  C  is  higher  than  that  of 
A,  and  the  decrease  in  the  tensile  elastic  modulus  is  very  slight, 
and  the  interface  strength  of  composite  materials  B  and  C  are 
higher  than  that  of  A,  therefore,  the  transverse  impact  strength 
of  composite  materials  B  and  C  is  higher  than  that  of  A. 

The  fracture  process  of  a  composite  material  longitudinal 
impact  specimen  is  relatively  complex.  Because  the  total  energy 
absorbed  when  the  specimen  is  broken  (not  including  the  energy 
consumed  as  the  specimen  flies  away)  includes  the  energy  consumed 
in  the  "initial  fracture"  U^,  and  the  fracture  work  consumed  in 
the  "crack  propagation"  process  Up,  therefore,  it  is  necessary  to 
comparatively  analyze  the  Ui  and  Up  of  the  three  composite 
materials  in  order  to  understand  the  contribution  of  the  matrix 
characteristics  to  the  impact  resistant  properties  of  composite 
materials . 
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Figure  4.  Failed  transverse  impact  specimen 

Because  the  interface  strength  of  the  three  composite 
materials  increases  with  increasing  matrix  toughness ,  therefore , 
the  capabilities  in  matrix  transfer  and  load  balancing  also  in-  / 
crease  with  increasing  toughness.  Furthermore,  the  tensile 
strength  and  fracture  elongation  ratio  of  the  toughened  matrices 
B  and  C  are  higher  than  those  of  matrix  A.  Consequently,  a  rel¬ 
atively  larger  stress  is  needed  for  the  destruction  of  the  micro¬ 
scopic  dimension  in  the  matrix  or  between  the  fibers  to  take 
place.  Therefore,  the  toughening  of  the  matrix  increases  the 
"initial  fracture"  work  Uj_  of  the  composite  material. 

An  observation  of  failed  longitudinal  impact  specimens  made 
of  composite  materials  A,  B,  and  C  discovered  that  the  mode  of 
fracture  of  these  three  materials  is  mainly  due  to  fiber  fracture. 
Only  a  small  amount  of  interface  detachment  existed  (see  Figure  5) 
This  is  because  all  the  epoxy-glass  fiber  composite  materials 
have  the  sufficient  interface  strength.  Hence,  it  is  not  pos¬ 
sible  to  create  a  great  deal  of  interface  detachment  in  the 


"fracture  propagation"  process.  Fiber  fracture  is  always  the 
major  mode  of  the  fracture  of  the  specimen  as  a  whole.  There¬ 
fore,  the  0  of  this  type  of  composite  material  must  increase 
P  rgi 

with  increasing  interface  strength.  J 

According  to  the  above  analysis,  the  presence  of  a  differ¬ 
ence  among  the  longitudinal  impact  strength  of  composite  mate¬ 
rials  A,  B,  and  C  may  be  primarily  due  to  the  difference  in  the 
toughness  of  the  matrix  materials  and  in  the  interface  strength. 
However,  because  there  are  more  factors  influencing  the  failure 
of  a  longitudinal  impact  specimen,  the  effect  of  the  matrix 
toughness  on  the  longitudinal  impact  strength  will  be  interfered 
with  by  other  factors. 

The  above  analysis  and  the  results  of  impact  strength  test¬ 
ing  showed  that  the  improvement  of  the  epoxy  resin  matrix  tough¬ 
ness  could  raise  the  resistance  against  fracture  of  a  composite 
material  in  a  high  speed  impact  state. 


Figure  5.  Failed  longitudinal  impact  specimen 
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Figure  6.  K-a  curves  of  unidirectional  E-glass  fiber  composites 

1.  kg.mm^ 

2 .  mm 

3 .  initial  crack  length  aQ=3 . 5mm 
(3)  Fracture  Toughness 

Similar  to  metallic  materials ,  the  stress  strength  factor 
K  can  be  used  as  a  fracture  mechanical  parameter  for  composite 
materials.  However,  different  from  metallic  materials,  it  is 


very  difficult  to  arbitrarily  choose  a  thickness  for  a  composite 
material  specimen.  Furthermore,  using  a  thick  specimen,  it  is 
very  difficult  to  ensure  that  a  crack  will  propagate  on  the  same 
plane.  It  is  not  easy  to  obtain  the  plane  strain  fracture  tough 
ness  corresponding  to  a  metallic  material.  Therefore,  in  this 
paper  the  stress  strength  factor  at  the  crack  instability  point 
of  a  specific  thickness  specimen  was  used  as  the  fracture  tough¬ 
ness  Kg.  The  instability  point  was  determined  by  the  K-a  curve 
for  crack  propagation.  The  K-a  curves  for  the  three  unidirec¬ 
tional  glass  fiber  composite  materials  A,  B,  and  C  are  shown  in 


Figure  6. 

Analysis  of  Figure  6:  the  vertical  lines  in  each  curve 
corresponds  to  the  first  stage  of  crack  propagation.  In  this 
stage,  the  crack  tip  displacement  increases  with  increasing  load. 
However,  the  crack  has  not  yet  begun  to  propagate.  The  crack 
length  is  a  fixed  value.  Immediately  afterwards,  it  enters  the 
second  stage,  a  sub-crack  region  begins  to  appear  at  the  tip  of 
the  crack.  With  an  increasing  stress  strength  factor,  the  crack 
strength  increases  slightly.  The  third  stage  of  crack  propaga¬ 
tion  is  the  gradual  stable  propagating  process  of  the  crack.  In 
this  stage ,  the  increment  of  crack  length  and  the  increment  of 
the  stress  strength  factor  showed  a  linear  relation.  The  first 
turning  point  on  the  K-a  curve  begins  to  appear.  When  the  stress 
strength  factor  increases  to  a  certain  value,  the  slope  of  the 
straight  line  section  changes.  The  second  turning  point  appears. 
This  is  the  critical  point  at  which  the  crack  changes  from  stable 
propagation  to  unstable  propagation.  The  straight  line  sections 
beyond  the  turning  point  points  upward  or  downward.  It  is  deter¬ 
mined  by  the  stress  transferred  by  the  fibers  from  the  crack  tip 
to  the  hardened  band. 

The  test  results  of  the  fracture  toughness  of  the  three  com¬ 
posite  materials  are  shown  in  Table  5.  In  the  paper,  Kq,  and 
are  the  stress  strength  factors  corresponding  to  the  crack  open¬ 
ing  point  and  the  instability  point,  respectively. 

From  the  data  shown  in  Table  5  one  can  see  that  the  value 
of  Kq  of  a  composite  material  increases  with  increasing  matrix 
toughness.  This  is  because  with  increasing  matrix  material 
toughness  the  damaged  region  at  the  tip  of  the  crack  created  by 
the  propagation  of  the  sub-critical  crack  becomes  larger.  Con¬ 
sequently,  more  energy  can  be  absorbed  and  the  concentration  of 
stress  is  reduced.  Therefore,  a  larger  stress  is  required  to 
create  a  crack.  The  data  in  Figure  5  also  shows  that  with  in¬ 
creasing  matrix  toughness, 
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Table  5, 


Fracture  toughness  of  unidirectional  S-glass  fiber 
composites 

1 .  epoxy  resin  matrix 

2.  composite  material 

3.  Kq,  Kg/mm3/2 

4.  Kq,  Kg/mm^/^ 
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Table  6.  Mechanical  properties  of  unidirectional  E-glass  fiber 
composites 

1 .  test  item 

2.  test  results  2 

3.  transverse  tensile  strength,  kg/cnr  5  2 

4.  transverse  tensile  elastic  modulus,  xlO3 kg/cni 

5.  longitudinal  tensile  strength,  kg/cm2  _  ~ 

6.  longitudinal  tensile  elastic  modulus,  xlCrkg/cni 

7.  static  flexual  strength,  kg/cnr  2 

8.  flexual  elastic  modulus,  xlO^kg/cnr  ~ 

9.  transverse  compression  strength,  kg/cni  2 

10.  longitudinal  compression  strength,  kg/cnr 

11.  short  beam  shear  strength,  kg/cm^ 


S3 


12.  glue  content,  % 

1 3 .  average 

14.  scatter  coefficient 

15.  average 

16.  scatter  coefficient 

17.  average 

18.  scatter  coefficient 

19.  average 

20.  scatter  coefficient 

21.  average 

22.  scatter  coefficient 

23.  average 

24.  scatter  coefficient 

25.  average 

26.  scatter  coefficient 

27.  average 

28.  scatter  coefficient 

29.  average 

30.  scatter  coefficient 

31.  composite  material 

the  value  of  of  the  composite  material  also  increases  corre-  / 
spondingly.  Composite  materials  C  and  B  are  8%  and  24%  higher 
than  A,  respectively.  However,  comparing  with  the  amplitude  of 
increase  in  the  toughness  of  the  matrix  material  C,  it  decreases. 
This  may  be  due  to  the  fact  that  the  matrix  material  in  a  fiber 
composite  material  is  under  the  action  of  forces  from  three  dir¬ 
ections.  The  tough  matrix  does  not  yield  easily,  and  consequently 
it  fails  by  embrittlement. 

2.  Conventional  Mechanical  Properties  of  Unidirectional 
Glass  Fiber  Composite  Materials 

The  major  mechanical  properties  of  composite  materials  A,  B, 
and  C  are  shown  in  Table  6 .  The  variation  of  each  property  with 
increased  toughness  is  shown  in  Figure  7. 

Figure  7  shows  that  after  the  matrix  is  toughened  the  trans¬ 
verse  tensile  strength,  transverse  compression  strength,  and 
short  beam  shear  strength  of  composite  materials  B  and  C  increase. 
The  transverse  tensile  modulus  decreases.  The  changes  of  other 
properties  are  not  obvious. 
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Figure  7.  Effect  of  matrix  toughness  on  the  mechanical  proper¬ 
ties  of  unidirectional  E-glass  fiber  composites. 

I.  Transverse  tensile  strength 

II.  Transverse  compressive  strength 

III.  Short  beam  shear  strength 

IV.  Transverse  tensile  modulus 

V.  Longitudinal  tensile  strength 

VI.  Longitudinal  tensile  modulus 

VII.  Flexural  strength 

VIII.  Flexural  modulus 

IX.  Longitudinal  compressive  strength 

1 .  relative  value 

Because  the  reinforcing  material  of  the  three  composite 
materials  is  the  same  type  of  glass  fiber,  therefore,  the  differ¬ 
ence  in  the  mechanical  properties  of  a  composite  material  is 
mainly  a  reflection  of  the  different  mechanical  properties  of 
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the  three  matrix  material  and  the  different  fiber-matrix  inter¬ 
face  strength. 

The  data  shown  in  Table  2  indicated  that  the  major  differ¬ 
ence  in  the  convention  mechanical  properties  of  the  matrix 
materials  A,  B,  and  C  is  shown  in  terms  of  tensile  strength, 
fracture  elongation  rate,  and  tensile  elastic  modulus.  The  first 
two  properties  increase  with  increasing  matrix  toughness.  The 
latter  decreases  slightly.  From  the  short  beam  shear  strength 
of  a  composite  material  one  knows  that  the  interface  strength  of 
composite  materials  B  and  C  is  higher  than  A.  Furthermore,  C  is 
higher  than  B. 

Because  of  the  aforementioned  difference  present  in  the 
tensile  properties  of  the  matrix  materials  A,  B,  and  C  and  inter¬ 
face  strength  of  a  composite  material,  therefore,  the  convention¬ 
al  mechanical  properties  of  the  three  composite  materials  also 
vary  correspondingly.  Also  because  the  toughness  of  the 
matrix  is  directly  related  to  the  tensile  properties  and  the  in¬ 
terface  strength  of  the  composite  material,  therefore,  the  mechan¬ 
ical  properties  of  the  composite  material  must  be  affected  by  the 
matrix  toughness.  However,  the  extent  of  effect  on  various 
mechanical  properties  is  not  the  same . 

IV.  CONCLUSIONS 

As  the  matrix  toughenss  increases,  the  fracture  resistant 
properties  (impact  strength,  fracture  surface  energy,  fracture 
toughness)  of  an  epoxy  resin-unidirectional  glass  fiber  compos¬ 
ite  material  also  increase  correspondingly.  The  short  beam 
shear  strength,  transverse  tensile  strength,  and  transverse  com¬ 
pressive  strength  of  a  composite  material  also  increase  with 
rising  matrix  toughness.  But,  the  transverse  tensile  modulus 
decreases.  The  effect  of  matrix  toughness  on  the  flexural  proper¬ 
ties  and  longitudinal  mechanical  properties  is  not  apparent. 
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ABSTRACT 

A  dynamic  analysis  monitoring  technique  is  presented  for  autoclave  cu¬ 
re  process  of  composites  constructed  of  T300/8DTT  carbon-fiber  epoxy  resin 
pre-preg.  On  the  basis  of  the  monitoring  technique  and  with  a  domestic 
monitor  model  JK2107,  the  appropriate  time  of  pressure  application  can  be 
chosen  and  then  the  ‘optimum  band  of  pressure  application*  can  be  deter¬ 
mined  successfully  by  a  small  number  of  experiments. 

During  cure  process  of  thermosetting  resin  system,  a  dielectric  loss  tan¬ 
gent  vs  time  curve  caa  be  obtained  with  the  monitor  and  the  ‘optimum 
band  of  pressure  application*  can  be  found  from  the  curve.  Thus  it  is  po¬ 
ssible  to  choose  the  proper  time  of  pressure  application  by  the  loss  tan?e- 
nt«  and  composite  plates  with  a  low  void  content  are  made  by  monitoring 
the  time  chosen. 

This  paper  proposes  a  new  way  to  move  the  ‘optimum  band  of  pressu¬ 
re  application’  to  the  peak  of  loss  tangent  curve  or  behind  it  by  incre¬ 
asing  frequency  given  by  monitor  JK2107  and  to  use  the  peak  of  the  curve 
as  a  reference  signaL  Therfore*  the  time  of  pressure  application  can  be 
controlled  exactly  and  the  quality  of  composite  plates  is  guaranteed  steadily 
in  this  way. 
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Monitoring  Cure  Process  of  Composites  with  Dynamic 
Dielectric  Analysis  Technique 

Tu  Shicun,  Li  Xingzhou,  Shen  Rongying,  and  Man  Shupeng 
(Beijing  Aeronautical  Material  Institute) 

INTRODUCTION 

A  composite  material  using  a  high  modulus  carbon  fiber  as 
the  reinforcing  body  has  a  much  higher  specific  strength, 
specific  modulus,  and  fatigue  resistance  as  compared  to  those 
of  steel,  aluminum,  and  titanium  alloys.  As  far  as  reducing  the 
structural  weight  of  an  airplane  is  concerned,  this  is  a  type  of 
aircraft  structural  material  with  a  bright  future  in  development. 
According  to  the  prediction  of  experts ^ ,  the  composite  material 
used  in  the  high  performance  fighter  structure  will  rise  to 
around  fifty  percent  by  the  end  of  this  century  from  a  few  per¬ 
cent  at  the  present  time.  However,  the  manufacturing  of  parts 
by  composite  materials  frequently  causes  the  properties  of  the 
materials  to  fluctuate  because  there  is  an  absence  of  a  rigorous 
quality  control  method.  The  scatter  coefficient  of  the  material 
strength  is  increased.  Thus,  it  is  necessary  to  increase  the 
safety  coefficient.  The  superiority  of  composite  materials 
cannot  be  fully  materialized. 

Void  is  a  frequently  observed  defect  in  fiber  reinforced 
composite  material.  In  the  meantime,  it  is  also  one  of  the 
major  reasons  causing  the  properties  of  a  material  to  scatter. 
This  defect  causes  the  inside  of  the  composite  material  to  be¬ 
come  porous.  According  to  theoretical  derivation,  as  well  as 
large  amount  of  experimental  data,  the  presence  of  voids  was 
proven.  It  causes  the  mechanical  properties  of  the  composite 
material  to  decline.  In  the  seventies,  the  technical  workers 
in  and  out  of  our  country  were  extremely  concerned  with  the 
effect  of  voids  on  fracture,  fatigue,  shear  strength,  and  hygro- 
scopicity.  A  great  deal  of  work  was  carried  out  on  non-destruc¬ 
tive  techniques  to  detect  voids.  A  lot  of  articles  have  been 
published^ 2-4 ^ .  As  far  as  the  dynamic  mechanical  properties 

are  concerned,  a  void  is  a  stress  concentrating  point.  Under 


125  a 


the  action  of  a  repeated  load,  it  is  a  fatigue  source.  As  far 
as  static  mechanical  properties  are  concerned,  voids  cause  the 
inside  of  the  material  to  become  porous.  They  also  cause  the 
mechanical  properties  to  decline.  According  to  a  report ^  , 
in  an  epoxy  resin  based  carbon  fiber  composite  material,  when 
the  porosity  is  under  4%,  the  laminar  shear  strength  of  the 
composite  material  will  decrease  by  4%  to  7%  as  the  porosity 
increases  by  1%.  With  respect  to  other  resins,  certain  key 
properties  also  decrease  to  the  same  extent.  Some  of  them  even 
decrease  by  larger  amounts.  Therefore,  experts  have  suggested ^ 
that:  as  major  structural  components  of  an  aircraft,  it  is  re¬ 
quired  that  pores  should  not  be  present  in  composite  material 
fabricated  parts ,  or  they  should  be  rigorously  controlled  to 
under  1%.  As  far  as  secondary  structural  parts  are  concerned, 
the  porosity  should  not  exceed  2%.  From  the  above,  one  can 
clearly  see  the  significance  of  obtaining  low  porosity  composite 
material  fabricated  parts. 

In  this  paper,  a  non-destructive  detection  technique  which 
is  capable  of  obtaining  low  porosity  composite  materials  -  the 
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dynamic  dielectric  monitoring  technique1  ,  is  introduced. 

This  technique  can  be  used  to  monitor  the  solidification  process 
of  carbon  fiber  composite  materials. 

EXPERIMENTAL 

1.  Brief  Description  of  the  Method. 

The  so-called  dynamic  dielectric  analysis  technique  is  a 
measuring  technique  which  is  capable  of  automatically  and  con¬ 
tinuously  monitoring  the  variation  of  dielectric  properties  of 
resin  systems  during  solidification  under  programmed  temperature 
control . 
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The  dielectric  analysis  technique  uses  the  motion  of  polar  /70 
groups  present  in  a  thermosetting  resin  to  express  the  state  of 
the  resin  system  at  a  particular  instance  in  the  solidification 
process.  When  an  alternating  and  varying  voltage  is  applied  to 


the  specimen,  due  to  the  interaction  between  the  dipole  moment 
and  the  electric  field,  the  dipole  moment  will  rotate  along  the 
direction  where  the  electric  field  is  generated.  Because  these 
polar  groups  are  attached  to  other  molecules,  therefore,  their 
movement  is  limited.  Hence,  the  extent  to  which  dipole  moment 
can  rotate  freely  is  related  to  the  extent  to  which  these  polar 
groups  can  move  freely  in  the  polymer.  During  a  polymerization 
reaction,  the  result  is  always  to  increase  the  constraints  ap¬ 
plied  to  the  migration  of  dipole  moments.  Therefore,  the  di¬ 
electric  response  characteristics  of  a  resin  system  with  respect 
to  an  applied  external  alternating  and  varying  electric  field 
will  vary  as  the  reaction  process  proceeds.  The  alignment  of  a 
dipole  moment  and  the  lagging  effect  when  it  rotates  can  be  ex¬ 
pressed  by  the  dielectric  coefficient  and  the  dielectric  loss 
factor,  respectively.  Furthermore,  both  of  them  are  functions 
of  frequency  and  temperature.  The  dielectric  monitoring  tech¬ 
nique  uses  the  variation  patterns  of  these  two  physical  quanti¬ 
ties  to  monitor  and  control  the  relevant  technological  parameters 
The  curves  of  dielectric  characteristics  (capacitance  C  and  the 
tangent  of  the  dielectric  loss  angle  tgfi )  and  viscosity  ( V  )  in 
a  typical  epoxy  resin  solidification  process  are  as  shown  in 
Figure  1.  From  Figure  1  one  can  find  the  portions  corresponding 
to  the  three  stages  in  the  solidification  process.  In  the  solid¬ 
ification  process,  the  temperature  (T)  and  the  polymerization 
reaction  affect  the  direction  of  the  dipole  moment  at  various 
degrees  alternatively.  During  the  initial  stage  of  solidifica¬ 
tion,  the  resin  begins  to  soften  with  increasing  temperature. 

It  is  gradually  changed  into  a  molten  state.  The  activity  of 
the  polar  group  increases  with  decreasing  viscosity.  The  elec¬ 
trical  capacitance  also  continues  to  increase.  Due  to  the  evap¬ 
oration  of  low  molecular  weight  compounds,  the  energy  loss  con¬ 
tinues  to  rise.  A  flow  peak  is  shown  in  the  loss  curve.  As 
the  viscosity  of  the  resin  further  decreases,  it  will  reach  its 
minimum . 


Figure  1.  Change  in  viscosity  and  dielectric  chaxacteristic 
during  curing  of  typical  epoxide  resin. 

This  corresponds  to  the  first  stage  in  the  solidification  process. 
Then,  the  polymerization  reaction  begins.  It  causes  the  viscos¬ 
ity  to  increase.  Moreover,  it  hinders  the  alignment  of  a  dipole 
moment  along  an  electric  field.  Consequently ,  it  causes  tgfi  to 
once  again  increase.  The  electrical  capacitance  also  decreases. 
Due  to  the  furthering  of  solidification,  the  cross-linking  re¬ 
action  causes  the  viscosity  to  increase  rapidly.  At  this 
moment,  tg 6  decreases  from  its  peak  and  gradually  levels  off. 

This  means  that  the  solidification  has  basically  been  completed. 
Therefore,  we  can  monitor  and  control  the  solidification  process 
on  a  real  time  basis  based  on  the  flow  variation  characteristics 
of  the  resin  by  using  the  dielectric  analysis  technique. 

2 .  Experimental  Apparatus 

In  this  experiment,  a  domestic  Model  JK  2107  automatic  di¬ 
electric  monitoring  device,  which  is  capable  of  automatically 
reading  and  continuously  recording  T,  C,  and  tgS  was  used. 

The  measuring  range  of  this  device  included:  0.1 ~ 10  KHz  in 
frequency,  0  ~1000  pF,  0~  «>  for  tgS  ,  and  0~  350°C  in  tempera¬ 
ture.  The  block  diagram  of  this  dielectric  monitoring  device 
is  shown  in  Figure  2. 


128 


Figure  2.  Monitoring  curve  of  dielectric  method. 

1.  The  Model  JK  2107  Device (3)  (the  temperature 
program  control  and  measuring  part) 

2.  The  Model  JK  2107  Device  (the  dielectric  measur¬ 
ing  part) 

3.  The  Model  LZ3-304  Function  Recording  Device 

In  the  figure,  the  detector  (1)  is  a  parallel  plate  capacitor 
using  a  carbon  fiber  pre-preg  material  as  the  medium.  The 
entire  detector  is  placed  in  a  heated  pressurized  container  (2) . 
When  the  testing  frequency  is  selected,  using  a  certain  temper¬ 
ature  program,  the  Model  JK  2107  dielectric  instrument  (3)  is 
used  to  measure  and  output  the  information  regarding  to  tgfi ,  C , 
and  T  in  the  solidification  process  of  the  resin  system.  These 
data  points  are  used  as  time  functions  of  a  three  channel  scan¬ 
ning  recorder.  The  LZ3-304  function  recorder  (4)  is  used  to 
plot  the  correlation  curves  between  tg5 ,  C,  and  T  with  time  t. 

3.  Selection  of  Optimized  Technological  Conditions  and 
the  Monitoring  and  Control  Method. 

For  the  usual  composite  material  manufactured  parts,  the 
frequently  used  solidification  process  obeys  a  standard  set  of 
temperature,  time  and  pressure.  We  call  this  standard  as  the 
static  technological  quality  control  method.  Under  the  situa¬ 
tion  that  various  factors  are  rigorously  controlled,  it  is 
possible  to  obtain  very  good  results.  In  reality,  in  the  middle 
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of  a  complicated  matter,  once  there  is  an  interference  from 
ternal  factor  the  properties  of  the  fabricated  parts  will  de¬ 
crease  due  to  improper  solidification.  They  may  even  have  to  be 
scrapped.  This  is  the  major  disadvantage  of  the  static  technol¬ 
ogical  quality  control  method.  Therefore,  the  development  of 
a  non-destructive  monitoring  and  control  technique  for  the  resin 
system  to  regulate  the  technological  process  on  a  real  time 
basis  becomes  very  important. 

The  dynamic  dielectric  analysis  technique  was  used  to  moni¬ 
tor  the  solidification  process  of  a  resin  system.  Furthermore, 
the  dielectric  characteristic  curves  related  to  the  resin  solid¬ 
ification  properties  were  used  as  the  reference  curves  in  the 
selection  of  optimum  technological  parameters.  Then,  the  curves 
were  used  as  a  basis  to  regulate  the  technological  process  on  a 
real  time  basis.  It  is  then  possible  to  obtain  the  optimum 
properties  attainable  by  the  material.  The  quality  of  the 
material  is  stable  and  reliable. 

In  order  to  obtain  the  optimum  technological  parameters, 
this  experiment  was  carried  out  using  the  following  method. 

Under  the  condition  that  the  temperature  program  was  given  by 
DSC,  the  pressure  application  time  was  the  key  problem.  For 
this  reason,  we  carried  out  the  solidification  test  in  a  heated 
pressurized  container.  The  speciment  used  a  pre-preg  formed  by 
8DT7  epoxy  resin  and  T300  carbon  fiber.  The  specimens  were  14 
layer  0  laminated  10  x  15cm  samples.  The  pressure  was  7Kg/cm  . 
An  instantaneous  pressure  application  method  was  used.  On  the 
tgfi  curve  in  Figure  1,  5  pressure  application  points  at  differ¬ 
ent  times.  A,  B,  C,  D,  and  E,  were  chosen.  The  solidified  com¬ 
posite  material  was  tested  in  order  to  determine  the  short  beam 
shear  strength  (SBSS) ,  bending  resistance  strength  (af) ,  and 
porosity  (V%).  Experimental  results  are  plotted  in  Figure  3. 


Figure  3.  Effect  of  different  pressure  application  time  on 
properties. 

1.  pressure  application  band 

According  to  these  experimental  results ,  we  can  obtain  the 
following  conclusions: 

(1)  The  pressure  application  time  cannot  be  too  early,  nor 
too  late.  If  the  pressure  is  applied  too  early  (point  A  in 
Figure  4) ,  the  resin  is  lost  excessively.  In  the  meantime,  due 
to  the  fact  that  low  molecular  weight  compounds  are  not  vola¬ 
tilized  sufficiently,  the  porosity  is  increased.  In  a  local 
area,  the  internal  voids  are  connected  together  to  lead  to 
local  delamination.  If  the  pressure  is  applied  too  late  (point 
E  in  Figure  4)  beyond  the  cure  point,  the  resin  is  already  cured 
Externally  applied  pressure  can  no  longer  press  the  resin  and 
the  fiber  to  a  solid  and  compact  entity.  Thus,  local  resin 
accumulations  will  be  formed  and  ompact  entity.  Thus,  local 
resin  accumulations  will  be  formed  inside  the  material.  Fur¬ 
thermore,  larger  voids  will  also  appear. 

(2)  When  the  pressure  is  applied  too  early  or  too  late,  it 
will  lead  to  increasing  porosity  in  the  composite  material. 

The  strength  decreases.  Therefore,  in  order  to  obtain  a  low 
porosity  composite  material,  it  is  necessary  to  fully  remove 
low  molecular  weight  volatile  compounds.  Furthermore,  pressure 
must  be  applied  at  a  suitable  time  (point  C  in  Figure  3)  before 


it  starts  to  cure  so  that  the  optimum  properties  of  the  materials 
can  be  attained. 

(3)  From  Figure  3  one  can  see  that  a  pressure  application 
belt  can  be  found  in  a  region  around  point  C  in  which  the 
strength  variation  is  not  too  much.  Its  width  can  be  determined 
by  the  allowed  fluctuation  range  of  the  material . 
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Point  C  of  prcssu*  application  Point  E  of  preoar*  application 


Figure  4.  Photo  micrograph  of  different  time  of  pressure  /72 

application. 

(4)  The  tangent  of  the  dielectric  loss  angle  tg<5  in  the 
region  corresponding  to  the  pressure  application  zone  is  our 
information  to  monitor  and  control  the  application  of  pressure. 

The  usual  value  of  tg<5  is  affected  by  many  factors.  To  simply 
use  tgS  to  monitor  and  control  the  pressure  application  time 
cannot  reach  good  results.  However,  the  peak  value  of  tg<5 
represents  a  special  turning  point. 

The  authors  presented  a  method  to  increase  the  testing  fre¬ 
quency  so  that  this  turning  point  is  displaced  to  within  or  in 
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front  of  the  monitoring  and  control  region  {See  Figure  5) .  It 
is  used  as  the  monitoring  and  control  signal.  Thus,  the  ac¬ 
curacy  of  monitoring  and  control  can  be  increased. 


Figure  5.  Monitoring  curve  of  dielectric  method. 
1.  pressure  application  band 


CONCLUSIONS 

We  used  the  Model  JK  2107  instrument  to  monitor  and  control 
the  pressure  application  of  a  carbon  fiber  composite  material 
during  the  solidification  process  with  a  press  in  a  heated  and 
pressurized  container.  The  test  has  already  proven  that:  this 
monitoring  and  control  technique ,  when  used  in  the  quality  con¬ 
trol  in  the  solidification  process,  can  obtain  a  composite 
material  plates  with  less  than  1%  porosity. 
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Brief  News  of  Technical  Meetings  of  the  Society 

•jfc  The  Society,  together  with  the  Astronautics  Society  of 
China  and  China  Aerodynamics  Research  Society,  held  the  "Tech¬ 
nical  Discussion  Meeting  on  Viscous  Flow  and  Surface  Layer"  on 
October  31st  through  November  6th  in  1982  in  the  city  of  Fuzhou. 

Sixty  four  representatives  attended  the  meeting.  The  meeting 
received  35  technical  papers,  and  22  were  read  in  the  meeting. 

The  Society  held  nine  technical  exchange  meetings  be¬ 
tween  November  and  December  1982.  They  included: 

1 .  On  November  2-8 ,  the  "Technical  Meeting  on  Environmen¬ 
tal  Control  and  Protection  in  Aircraft  Cockpits"  was  held  in 
Kunming.  92  representatives  attended  the  meeting.  65  technical 
papers  were  received.  50  of  them  were  read  in  the  meeting. 

2.  On  November  5-11,  the  "Forging  Technique  Exchange  Meet¬ 
ing"  was  held  in  Jiujiang.  59  representatives  attended  the 
meeting.  42  technical  papers  were  received.  28  of  them  were 
read  in  the  meeting. 

3.  On  November  14-20,  the  "Technical  Exchange  Meeting  on 
the  Auxiliary  Engine  and  Precision  Machining  Special  Fields" 
was  held  in  Wuxi,  Jiangsu.  122  delegates  attended  the  meeting. 

It  received  85  technical  papers  and  57  of  them  were  read  in  the 
meeting . 

4.  On  November  15-21,  the  "Technical  Discussion  Meeting  on 
the  Strength  of  Aeronautical  Gears"  was  held  in  Zhuzhou,  Hunan. 

73  delegates  attended  the  meeting.  45  technical  papers  were 
received  and  34  of  them  were  read  in  the  meeting. 

5.  On  November  18-22,  the  "Rocket  Engine  Special  Field 
Technical  Exchange  Meeting"  was  held  in  Changsha.  47  deleaates  attended 
the  meeting.  40  technical  papers  were  received,  and  25  were  read. 

6.  On  November  22-28,  the  "Aircraft  Incident  Inspection 
and  Analysis  Technical  Discussion  Meeting"  was  held  in  Beijing. 

73  delegates  attended  the  meeting.  35  technical  papers  were 
received,  and  25  of  them  were  read  in  the  meeting. 

7.  On  December  4-10,  the  "Aeronautical  Auxiliary  Engine 
Computer  Aided  Design  and  Manufacture  Meeting"  was  held  in 
Fuzhou.  79  representatives  attended  the  meeting.  41  technical 


papers  were  received,  and  36  of  them  were  read  in  the  meeting. 

8.  On  December  6-10,  the  "Technical  Exchange  Meeting  on 
Data  Collection,  Processing,  and  Control  Systems  Centered 
Around  Computers"  was  held  in  Chengdu.  79  delegates  attended 
the  meeting.  50  technical  papers  were  received,  and  46  of  them 
were  read  in  the  meeting. 

9.  On  December  16-23,  the  "High  Temperature  Fatigue  and 
High  Temperature  Strain  Testing  Technique  Meeting"  was  held  in 
Chengdu.  102  delegates  attended  the  meeting.  60  technical 
papers  were  received,  and  48  of  them  were  read  in  the  meeting 

The  Society  of  Aeronautics  in  China  Held  an  Academic 
Activity  Symposium  in  Beijing. 

The  Society  of  Aeronautics  in  China  held  an  academic 
activity  symposium  on  January  18,  1983  at  the  People's  Confer¬ 
ence  Hall  in  Beijing.  There  were  ninety  one  people,  including 
the  secretary  general  and  vice  secretary  general  of  the  Society 
of  Aeronautics  in  China  in  Beijing  together  with  the  chairmen 
of  various  special  field  committees,  academic  secretaries,  and 
responsible  comrades  of  related  organizations .  The  secretary 
general  of  the  Society  of  Aeronautics  in  China,  comrade  Wang 
Junkui,  presided  the  meeting.  Comrade  Zhao  Zhengyan,  the  vice 
secretary  general  of  the  Society  of  Aeronautics  in  China ,  gave 
a  working  report  on  the  "summary  of  academic  activities  in  1982 
and  the  "key  points  of  academic  activities  in  1983."  Comrade 
Liu  Tianfu,  who  is  the  vice  secretary  of  the  Party  Committee 
of  the  Chinese  Science  Association,  arrived  at  the  meeting  and 
gave  an  important  speech. 

Under  the  leadership  of  the  Party  organizations  of  the 
Chinese  Science  Association,  Aeronautical  Industry  Ministry, 
and  the  Society,  and  with  the  support  of  the  vast  scientific 
and  technical  personnel ,  the  related  units  nationwide ,  and  the 
aeronautic  scoieties  in  various  provinces  and  cities,  through 
the  joint  effort  of  various  special  committees  and  special 
academic  groups,  the  academic  activities  of  the  Society  in  1982 
were  successfully  completed  according  to  the  1982  academic 
activity  plan.  Twenty  eight  academic  meetings  were 
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held  in  1982.  A  short  term  class  was  held.  As  compared  to  last 
year,  the  number  of  activities  increased  by  26%.  The  total  num¬ 
ber  of  technical  papers  received  was  1699,  and  1259  of  them  were 
read  in  various  meetings.  2499  people  participated  in  the  acad¬ 
emic  activities.  Through  these  academic  activities,  the  results 
of  recent  academic  research  and  working  experience  were  summar¬ 
ized  and  exchanged.  The  trends  of  development  of  related  special 
fields  in  and  out  of  our  country  were  analyzed  and  explored. 
Technical  advisory  suggestions  were  provided  to  the  relevant 
department.  The  various  academic  organizations  were  established 
and  further  strengthened.  The  future  academic  activity  plan  was 
discussed.  Through  these  academic  activities,  the  planning  of 
domestic  academic  activities  and  the  organizational  management 
work  were  further  strengthened.  The  rules  and  regulations  need¬ 
ed  were  established  and  perfected.  The  "temporary  regulation 
regarding  organizing  domestic  academic  meetings"  and  the  "brief 
rules  governing  special  committees."  The  work  of  various  spec¬ 
ial  committees  and  the  development  of  academic  activities  were 
promoted . 

The  1983  academic  activity  plan  has  already  been  formulated. 
A  total  of  35  academic  meetings,  5  training  classes,  and  2  acad¬ 
emic  working  meetings  will  be  arranged.  The  total  number  is  42. 
In  the  activities  in  this  year,  the  emphasis  will  be  placed  on 
the  selection  and  recommendation  of  outstanding  papers.  Condi¬ 
tions  will  be  created  so  that  training  classes  for  scientific 
and  technical  personnel  can  be  operated  well.  The  special  fea¬ 
tures  of  academic  activities  will  be  combined  in  order  to  open 
up  various  technical  consulting  activities.  The  extension  and 
exchange  of  academic  activity  results  must  be  seriously  carried 
out. 

Vice  Secretaries  of  the  Society,  comrades  Dong  Xiao,  Hang 
Mingyang,  and  Zhang  Peijian,  as  well  as  Special  Committee  Chair¬ 
men  such  as  comrades  Zhuang  Fengan,  and  Lin  Shier,  spoke  at 
the  meeting.  They  provided  valuable  suggestions  regarding  the 
direction  of  future  development  of  the  academic  activities  of 
the  Society  and  the  development  of  the  aeronautic  industry  in 


our  country.  During  the  meeting,  the  experts  and  professors  in 
the  elder  generation  of  the  aeronautic  community  in  our  country 
and  the  young  and  middle  age  technical  workers  gathered  together 
in  the  same  room  in  harmony.  In  the  afternoon,  the  various 
special  committee  academic  secretaries  in  Beijing  continued  to 
meet  specifically  with  respect  to  some  actual  problems  in  acad¬ 
emic  work.  Everyone  was  enthusiastically  in  agreement  with  the 
academic  working  conference  held  by  the  Society  in  August  this 
year.  The  academic  work  experience  was  summarized  and  exchanged. 
Problems  such  as  how  to  further  improve  the  quality  of  academic 
meeting  was  explored. 

Obituary  In  Deep  Memory  of  Comrade  Wang  Terong,  /77 

Chief  Editor  of  {(Journal  of  Aeronautics)) 

Chief  Editor  of  ( (Journal  of  Aeronautics) ) ,  Comrade  Wang 
Terong  died  on  December  2,  1982  due  to  sickness.  He  was  seventy 
four.  His  departure  was  an  important  loss  to  the  research  and 
education  in  aeronautics  in  our  country.  It  is  also  an  import¬ 
ant  loss  to  this  publication. 

Comrade  Wang  Terong  graduated  from  Shanghai  Jiaotong  Uni¬ 
versity  in  1932.  He  went  to  England  to  study  in  1935.  He  re¬ 
turned  to  our  country  in  1938.  Until  1951,  he  was  professor 
and  chairman  at  Southwest  Consolidated  University  and  Qing  Hua 
University.  In  1952,  he  participated  in  the  planning  work  of 
Beijing  Institute  of  Aeronautics  and  Astronautics.  For  over 
thirty  years,  he  had  made  great  contributions  to  the  construc¬ 
tion  and  growth  of  Beijing  Institute  of  Aeronautics  and  Astro¬ 
nautics.  In  1964,  during  the  founding  process  of  the  Society  of 
Aeronautics  in  China,  he  was  one  of  the  active  supporters.  Since 
1965,  from  the  beginning  of  ((Journal  of  Aeronautics)),  he  was 
the  chief  editor.  He  worked  very  hard  in  the  manuscript  review¬ 
ing,  editing,  and  publishing  of  each  volume  of  the  ((Journal)). 

Even  when  he  was  critically  ill,  he  still  remembered  the  work  at 
the  ((Journal)). 

For  over  half  a  century,  comrade  Wang  Terong  dedicated  all 
his  energy  to  the  aeronautic  industry  in  his  own  country. 


During  the  period  ruled  by  the  reactionary  Nationalists,  he  re¬ 
mained  straight  while  facing  the  dark  ruling  power.  He  did  not 
try  to  get  close  to  the  center  of  power.  He  was  dissatisfied 
with  the  corrupt  authority.  Furthermore,  he  was  actively  helping 
progressive  students  with  sympathy.  In  1949,  Comrade  Wang  Terong 
was  touring  the  United  States .  As  soon  as  he  heard  the  good 
news  that  the  new  China  was  founded,  he  immediately  started  to 
return  to  our  country  to  devote  to  the  construction  of  our  home¬ 
land.  He  was  actively  close  to  the  party  organization  and  in¬ 
sisted  on  the  accurate  position.  He  was  persistent  and  maintain¬ 
ed  the  realistic  attitude.  He  supported  the  leadership  of  the 
Party  and  the  socialistic  route.  During  the  three  year  diffi¬ 
culty  period,  he  not  only  did  not  withdraw,  on  the  contrary,  ob¬ 
tained  more  progress  politically  with  enthusiasm.  In  December 
1961,  he  joined  the  Chinese  Communist  Party  with  honor.  During 
the  ten  year  period  of  chaos,  no  matter  whether  it  was  political 
oppression  or  torture  in  life,  the  Comrade  Wang  Terong' s  decision 
to  struggle  for  the  Party  and  the  people  was  not  shaken.  After 
Premier  Chou  proposed  the  magnificant  objective  of  the  Four 
Modernizations,  he  gave  lectures  to  research  organizations  ac¬ 
tively.  He  provided  a  lot  of  rationalized  opinions  to  overcome 
technical  obstacles.  He  solved  some  major  difficulties  in  re¬ 
search  and  production.  Until  he  was  critically  ill,  he  still 
made  his  best  effort  to  illuminate  the  scientific  exploration 
route  with  his  life. 

Throughout  his  life.  Comrade  Wang  Terong  devoted  his  honor¬ 
able  life  to  the  aeronautical  industry  in  our  country.  It  was 
a  life  seeking  for  progress.  It  was  a  life  of  hard  struggle 
and  hard  work.  He  had  raised  many  outstanding  students.  He  was 
an  elder  leader  and  an  elder  teacher  with  reputation  in  the 
aeronautical  community.  However,  he  never  sought  for  honor  and 
position  for  himself.  He  was  a  famous  professor  and  an  expert. 

He  was  a  council  member  of  the  Chinese  Mechanics  Society,  a 
council  member  of  the  Society  of  Aeronautics  in  China,  a  speci¬ 
ally  invited  committee  member  of  the  Scientific  Technical  Com¬ 
mittee  of  the  Aeronautical  Industry  Ministry,  and  a  member  of 


■V- 


r 


the  Political  Association  of  Beijing.  However,  he  never 
claimed  credit  and  became  arrogant.  He  was  humble  and  easy  to 
approach.  He  had  written  and  reviewed  a  great  deal  of  papers. 
However ,  he  never  let  each  detail  go  and  repeatedly  pondered  on 
the  key  words.  Comrade  Wang  Terong  was  indeed  an  outstanding 
teacher  of  the  people.  He  was  indeed  a  good  member  of  the 
Chinese  Communist  Party. 

Comrade  Wang  Terong  left  us.  This  made  us  sincerely  sorrow¬ 
ful.  However,  we  must  transform  sorrow  into  strength.  We  must 
learn  from  him  the  excellent  quality  of  his  loyalty  to  the  Party 
and  to  education.  We  must  learn  the  rigorous  and  serious  atti¬ 
tude  in  studying  from  him.  Under  the  guidance  of  the  spirit  of 
the  12th  Congress  of  the  Party,  we  must  succeed  the  career 
started  by  Comrade  Wang  Terong.  We  must  continue  to  work  hard 
to  run  ((Journal  of  Aeronautics))  well  in  order  to  struggle  hard 
for  the  beautiful  tomorrow  of  the  aeronautical  and  astronautical 
career  of  our  country. 
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Editing  Committee  of  ( (Journal  of 
Aeronautics) ) 
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